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The Effect of Microgravity on the Growth of Lead Tin Telluride 
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Summary: The main objective of this research was to present a model for the 
prediction of the effect of the microgravity environment on the growth of 
Lead Tin Telluride. 

The attitude change and its relation to the experimental objectives: The 

main objective for the AADSF experiment on USMP 3 involving LTT growth 
was to estimate the effect of ampoule orientation on the axial and radial 
segregation of tin telluride. As the furnace was not situated on a gimbal 
there was no possibility to reorient the ampoule during the flight. Instead the 
only way to change the growth orientation was to change the attitude of the 
orbiter. This was accomplished by vernier rocket firings. 

In what follows it must be noted that the orbiter body coordinates are 
such that the positive z axis points outward from the belly', the positive x 
axis points outwards from the nose and the positive 'y' axis points outwards 
from the starboard side. The furnace which was in the pay load had its axis 
aligned with the orbiter's ' z' axis with the hot end closest to the shuttle 
body. There were basically three orientations that were desired. These 
corresponded to the ampoule being seen as a heated from above ( thermally 
stable-solutally unstable) configuration, the heated from below ( where the 
instabilities were reversed from the first orientation) configuration and an ' in 
between' case where the ampoule axis was misaligned with respect to the 
orbiters ' gz ' axis. 

In order to understand the role of each orbiter attitude that was 
requested it is necessary to have an idea of the first order effects on the 
residual acceleration levels on the orbiter. The center of gravity of the orbiter 
(CG) generally did not change very much during the USMP 3 part of the 
mission as the change in mass was insignificant. A displacement from the 
CG of the AADSF furnace and the attitude of the orbiter affected the so 
called “gravity gradient' force ( arising from the centrifugal component of 
acceleration). This gravity gradient had a major effect on the force fields 
which at low frequency levels ( under 1 Hz) was of the order of 10® ge 
where ge is earth's gravity. A second major effect on the residual forces is 
the atmospheric drag on the orbiter. This drag can only have a deceleration 
effect on the vehicle and therefore affect the forces on the particles in the 
AADSF. This drag was itself affected by the orbiter attitude, the position of 
the orbiter with respect to the sun { i.e., atmospheric density) and the 
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orbiting path around the globe as the earth is not a perfect sphere. The 
deceleration and the tendency for the orbiter to get into an aerodynamically 
stable mode necessitated the need for vernier booster firings. The forces 
associated with these corrective measures were of high magnitude ( roughly 
10 ^ge -10 ‘’ge) and were of a high frequency ( 5-10 Hz). 

Keeping in mind the main effects on the residual low frequency 
acceleration we now turn to the various attitudes that were chosen during 
the mission. The first attitude had a pitch of 185 degrees and a roll of 7 
degrees. This attitude is roughly equal to the orbiter flying in a position of 
payload to earth and tail into the wind. ( also called "Zlv ,"Xvv ). The drag was 
very small and the small roll angle served mainly to offer a minor change in 
the gravity gradient. The cold end of the furnace was closest to earth and 
the effect was to have a large axial to normal ratio of the axial vector. The 
second attitude was meant to give a heated from below configuration. At 
first sight it might appear that the best attitude would be a zero degree pitch. 
However this was considered to be a risky attitude on account of debris. 
The second attitude was therefore obtained by considering a deceleration 
mode of the orbiter. Here the orbiter was placed in a pitch of 90 degrees 
with a small negative yaw of 1 7 degrees ( in alternate terminology this 
amounted to an attitude of -Xiv and -i-Zw). The pitch virtually maximized the 
drag on the orbiter while the yaw served the purpose of adjusting the gravity 
gradient on the furnace. Because the orbiter was decelerating it had the 
opposite effect on the fluid particles in the AADSF ampoule relative to the 
shuttle and the net effect was to accelerate the fluid in the cold end towards 
the hot end. This was equivalent of generating a thermally unstable and 
solutally stable configuration. Because the deceleration force was not of a 
high magnitude the net axial to normal ratio was not very high. Moreover the 
normal component of the force fields were very small and so in the actual 
flight one could expect that small changes in the 'dead band' could cause a 
substantial undesired change in the axial to normal ratio. The last attitude 
had a pitch of 1 23 degrees and the net effect of the drag was to change the 
orientation of the axial vector in the furnace so that it behaved as if it was 
tilted with respect to gz of the shuttle. While the choice of orbiter attitudes 
had to do with the prefered growth direction the reasons for choosing these 
various growth directions or ampoule orientations resulted from preliminary 
work done using computational fluid dynamics. 

Preliminary estimates of axis orientation using CFD: The 

computational fluid dynamic calculation procedure is best explained by 
considering Figure 1 . The ampoule liquid region is assumed to be constant 
as the solidification rates are normally very small. The boundary conditions 
indicate hot and cold zones as well as insulating zones. The far field 
concentration is assumed to be constant and the interface condition respects 
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mass conservation along with solutal segregation. The effect of the 
magnitude of the gravity vector is seen in Figures 2 a) and b) and we 
immediately conclude that the flow at low frequency low amplitude 
accelerations (lO^ge) will be of a weak torroidal type. Higher amplitude 
forces will cause solutal convection to come into play but such high 
amplitude acceleration vectors were not present at the low frequency levels 
during USMP 3 and were mainly associated with high frequency activities 
such as booster firings and water dumps. While we do not presently have 
any predictions on the effects of high frequency g -jitter it is clear from the 
order of magnitude of the calculated velocity vectors for the case of low 
frequency with 1 0 ^ge that very little mixing takes place. The velocity is no 
greater than 10^ cm/ sec. Given an initial liquid region size of 5 cm, this 
small velocity amounts to an initial mixing time of 10® seconds. Meanwhile 
the solidification is at the rate of 1 cm/ hour. When the liquid region size is 

about 1 cm the mixing time is about 2 xIO ® seconds. Clearly this is 

insignificant because the entire growth period is about 2 xIO'* seconds . In 
other words we predict that only diffusion controlled growth ought to prevail 
at 10®ge and this more true at lower gravitational levels which were 
experienced during USMP 3. 

The effect of a five degree offset with respect to the vertical 
orientation was calculated and the results are graphically shown in Fig. 2 c. 
What is seen from this figure is that small tilts give rise to swirling flow and 
this flow contains the solutal boundary layer to the depleting surface. This 
may be contrasted with torroidal flow in Fig. 2 a ( for the vertical orientation) 
that sweeps the solute out of the solutal boundary layer. The solutal 

boundary layer contains most of the rejected SnTe and so swirling flow if 

anything should help by making diffusive growth more probable. In other 
words one might conclude that a constant off axis arrangement is better than 
an on axis aligned ampoule. Fig. 3 is a depiction of the mixing patterns that 
are seen when the ampoule is subjected to a time dependent tilt. It must be 
noted that the velocities are still very small and so even in the case when the 
tilt is a periodic function of time the growth is expected to be diffusion 
controlled. 

In summary we have concluded that diffusive growth was predicted 
under low frequency g- jitter conditions. The high frequency was not studied 
but we did conclude that the time constant for the fastest transporting 
mechanism ( heat transfer) was much larger than the corresponding period 
for high frequency ( 5 Hz) g- jitter. 

Publications: The publications that arose from this work and directly related 
to the objectives are attached and are listed below: 
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microgravity effects on convection and interfacial driven convection. These 
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ABSTRACT: 

A 3D numerical calculation is performed on a model that depicts the flow profiles due to thermo-solutal 
convection in a cylindrical tube. The calculations were done with the purpose of delineating the qualitative 
features of the flow profiles for the cases when the container's axis is perfectly aligned with respect to the 
mean gravity vector and also when it changes periodically with respect to the gravity vector. It is found that 
the flow profiles are similar to those of the Rayleigh-Benard problem in the case of perfect alignment while a 
swirling pattern appears when the tube’s axis is not aligned with the gravity vector. This indicates that it 
might be preferable to have a slight tilt in the container axis during crystal growth as svraiing flow will 
diminish axial mixing. The solutal convection is the dominant feature of the flow and is affected 
considerably by the gravity level. © 1999 COSPAR. Published by Elsevier Science Lid. 


INTRODUCTION 

This is a brief report describing the flow profiles that are induced in a low gravity environment in a 
Bridgman tube in which the fluid occupies a constant volume. The Bridgman tube as considered in this 
study is merely a circular cylinder that is subjected to radial thermal gradients and axial solutal gradients. 
Typically, the Bridgman tube is used in the vertical directional solidification of compound semi conductors 
such as Lead Tin Telluride. The growth of such materials is affected substantially by the convective flow 
profiles that accompany the process. This convection is due to thermal and solutal gradients that are 
generated because of the solidification process. Arnold eL al. (1991) did calculations to model a GaAs space 
experiment and concluded that three-dimensional flows occur under certain gravitational values and 
orientations. Their calculations were not concerned with solutal convection. Naumann and Baugher ( 1 992) 
have made analytical estimates of radial segregation in Bridgman growth for low-level steady and periodic 
accelerations. In any actual growth process, the liquid zone is ever shrinking and this can be expected to 
change the flow profiles quite a bit. Nonetheless it would be interesting to have an idea of the flow profiles 
that are generated when the force conditions on the ampoule are compatible to a time dependent 
microgravity level and where certain assumptions such as a constant liquid zone is assumed. We present here 
a numerical model that shows the effects of off axis and a time dependent orientation on the flow profiles in 
a Bridgman tube. The effect of tilting the otherwise vertical container with respect to gravity is also 
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described. The calculations were done with the purpose of delineating the qualitative features of the flow 
profiles for the cases when the container's axis is perfectly aligned with respect to the mean gravity vector 
and also when its axis periodically changes with respect to the gravity vector. The gravitational levels that 
are assumed range from lO'^ge (or ten micro g) to 10 ^gc where gc is earth's gravity. A value of ten micro g is 
reasonable as it is a fair representation of the low frequency accelerations experienced on the cargo bay of the 
U.S. space shuttle or on the future international space station if the Bridgman tube. A level of 100 micro g is 
not very probable, however we also present calculations that include this extreme case. To the best of our 
knowledge this is the only study that shows the dominant effect of solutal convection over thermal 
convection at the higher gravity levels and also the only study where the ampoule axis orientation is varied 
with time. 


Table 1 The Thermophysiceil Properties Used in the Calculations 


Density 

7.04 g/cm^ 

Kinematic Viscosity 

0.0024 cm^/s 

Thermal Diffusivity 

0.03 cm Is 

Solutal Diffusivity 

1 1 0'^ cm^/s 

Thermal Expansion Coefficient 

1.18 10"* /C° 

Solutal Expansion Coefficient 

0.22 /Weight fraction 

Segregation Coefficient 

0.7 


THE MODEL AND THE NUMERICAL SCHEME 

The model that is used assumes that the Boussinesq equations hold. Further the c 2 ilculations were done 
assuming that the fluid is Lead Tin Telluride reflecting our interest in compound semiconductors. The 
thermophysical properties of Lead Tin Telluride as used in the calculation are given in Table 1. Figure 1 
describes the situation when a container is subject to thermal gradients with a solidifying interface at 
z=H-Hsoiid. The thermal and concentration boundary conditions imposed on the container are given in Figure 
1. No-slip conditions are used at ail boundaries including the solid-liquid interface upon which the 
coordinate system is fixed and which is assumed to move down in the z direction at a constant speed Vs 
equal to 1 cm/hr. The height of the tube of diameter equal to 1 . cm.is given by H and assumed to be 5.0 cm., 
equally divided between the solid and liquid zones while the insulation zone is assumed to occupy the 
middle one third. These correspond roughly to the experimental ampoule used in a Lead Tin Telluride 
experiment that was conducted on USMP 3. The hottest temperature is assumed to be 1150 degrees Celsius 
while the coldest temperature is fixed to be 700 degrees, the interface being at 900 degrees. The orientation 
of the gravitational acceleration is expressed in terms of the angle between the gravity and the negative 
direction of the z-axis. The major assumption is that the liquid length is kept constant. As a result it is 
assumed that the end of the liquid was at a constant concentration , Co equal to 0.2 weight fraction. This is 
tantamount to a continuous feeding of such liquid at the solidification rate Vs . Before we go on it might be 
useful to point out that the thermal Rayleigh number is estimated to be about 65 for a gravity level lO’^ge 
while the solutal Rayleigh using the same length scale is about 147000 

















3D Numerical Model for Flow in a Bridgman Tube 


1399 



Fig. 1 The schematic of the geometry and thermal and 
concentration boundary conditions ofthe calculation. 


RESULTS OF THE CALCULATIONS AND CONCLUSIONS 

The finite volume method, SIMPLE (cf. Patankar,I980') was chosen to solve the governing equations. 
Figure 2 shows the flow profiles at two different gravity levels. Figure 2a describes the pattern that is 
expected at a constant g level of lO'^ge . What is to be observed is vertical stacking of an axisymmetric or 
torroidal pattern. This vertical stacking may be expected as the top of the ampoule is hotter than the bottom 
and the lower cell' is in the insulation zone. The configuration acts like a fluid that is heated from above' 
and the weak flow is primarily driven by radial gradients. The weakness in the lower cell is primarily due to 
the effect of the presence of the no slip' solid boundary. 

The situation changes somewhat for the case of a g level of lO'^ge as seen in Figure 2b, for here the solutal 
convection begins to play a part. The solutal gradients are unstable in the sense that they promote flow even 
if the thermal expansion coefficient is negligible. As observed earlier the solutal Rayleigh number is about 
147000 whereas the critical solutal Rayleigh number, in the absence of thermal gradients for this aspect ratio 
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(radius/ height =0.2) turns out to be about 50000 from the calculations of Hardin et. al.. The unicellular 
patterns that are observed are a result of the solutally driven convection and are expected for this geometry 
according to the results of the Rayleigh-Benard problem (cf, Hardin et.al., 1990). Notice further that the 
vertical stacking arrangement disappeared for the larger level indicating the dominance of the solutal 
convection over the thermal convection. Figure 2c shows the flow profiles resulting from a time invariant 
change of the ampoule axis with respect to the mean gravity vector of 10’ gc. What is immediately apparent 
is the swirling flow that helps contain the rejected solute near the solid liquid interface. The 'z' component of 
this type of flow is much smaller than the other wo components except near the interface and the end of the 
liquid where all components are set to zero on account of no-slip. It is a concentration induced flow because 
the velocity components near the interface in Figure 2c are much larger than in Figure 2a. In a real crystal 
growth configuration such an off axis tilt would help prevent axial mixing and therefore be beneficial to the 
crystal. If the 'g'- value is increased by an order of magnitude the flow is mostly of a unicellular style except 
near the interface. This is shown in Figure 2d. Figure 3 depicts the periodic change of the cylinder axis with 
the gravity vector. The frequency was set to be one cycle per hour. This was an arbitrary choice even though 
aerodynamic drag causes a readjustment every 20 minutes or so in a typical space orbiter. All the same the 
results can be expected to be qualitatively similar to those reported here in the case when the frequency is 
increased three fold. The flow profiles at a constant 'z' plane near the solid liquid interface are given at every 
quarter cycle, once a periodic steady state is reached. Notice that the direction of the swirl changes every 
half cycle i.e., when the gravity vector crosses the cylinder axis leading to local mixing near the solid liquid 
interface. 

The effect of gravity level and time periodic off axis alignment show collectively that the convection at low 
gravity in a bottom or top seeded Bridgman tube is primarily in the solutal driven mode as long as the gravity 
level is not very small. Moreover a slight tilt with respect to gravity causes the fluid flow to go into a 
swirling mode so that solute is contained near the solute-generating boundary. 
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Abstract 

One of the primary benefits of conducting scientific research in space is to take advantage of the low acceleration 
environment. For experimenters conducting space research in the field of materials science the quality of the science 
return is contingent upon the extremely low frequency acceleration environment ( < 1 Hz) aboard the spacecraft. Primary 
contributors to this low frequency acceleration environment (commonly referred to as the steady-state acceleration 
environment) include aerodynamic drag, gravity-gradient, and rotational effects. The space shuttle was used on the 
STS-75 mission as a microgravity platform for conducting a material science experiment in which a lead tin telluride alloy 
was melted and regrown in the Advanced Automated Directional Solidification Furnace under different steady-state 
acceleration environment conditions by placing the shuttle in particular fixed orientations during sample processing. The 
two different shuttle orientations employed during sample processing were a bay to Earth, tail into the velocity vector 
shuttle orientation and a tail to Earth, belly into the velocity vector shuttle orientation. Scientists have shown, through 
modeling techniques, the effects of various residual acceleration vector orientations to the micro-buoyant flows during 
the growth of compound semiconductors. The signatures imposed by these temporally dependent flows are manifested in 
the axial and radial segregation or composition along the crystal. 
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1. Introduction 

1,1. Motivation for conducting specific materials 
science experiments in space 

Crystal growth research in ground-based labo- 
ratories is complicated by the ever present effects of 
gravity acting on the sample during processing. 


* Corresponding author. 


Such an effect oftentimes drives natural convection 
within the sample. In theory, if these same types of 
experiments are conducted in a microgravity envi- 
ronment, buoyancy-driven convection effects due 
to gravity are essentially eliminated, thus providing 
ideal conditions for diffusive growth within the 
sample. Under a microgravity environment scien- 
tists are able to study the low frequency, low accel- 
eration (referred to as the steady-state acceleration 
regime) effects (if any) on materials. Furthermore, 
scientists can use a space platform to determine 
whether the particular direction of the steady-state 
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acceleration (SSA) vector present at the growth 
location influences the quality of the space-pro- 
cessed crystal. Results from these types of experi- 
ments will allow scientists to validate or update 
their models and theories in the field of crystal 
growth research. Such discoveries are expected to 
improve ground-based processing techniques for 
crystal processing [1-3]. 


2. Lead tin telluride grown in space 

2. 1. Description 

The lead tin telluride ( PbSnTe) alloy was selected 
as an experiment within the Advanced Automated 
Directional Solidification Furnace (AADSF) to fly 
aboard the Space Shuttle Columbia as part of 
the Third United States Microgravity Payload 
(USMP-3) mission (STS-75) in February and 
March of 1996. The material was selected as an 
experiment for investigation under microgravity 
conditions due to the fact that it can be used to 
produce infra-red detectors and lasers on Earth. 
The primary objective of the experiment was to 
determine the effects (if any) of the SSA environ- 
ment during the processing of each of the PbSnTe 
samples. Understanding the microscopic effects 
during processing in a low-g environment could 
improve the way these devices are made [4], 

2.2. Specific experimental studies of PbSnTe 

Even at acceleration levels as small as 5 x 10“ 
scientists interested in studying PbSnTe have de- 
vised theories on how the crystal microstructure 
will react to the small but constant accelerations 
present at the sample location [3]. Contributors to 
the SSA environment will be discussed later in the 
paper. 

Two possible crystal growth orientations with 
respect to the SSA (gravity) vector will be ad- 
dressed. Figs. 1 and 2 summarize these two possible 
orientations. Fig. I is described as the "hot on top" 
configuration; that is, the SSA vector travels in the 
direction from the hot temperature portion of the 
sample to the cooler temperature portion of the 


ig 


arripouie 

nnoves downwards 

Fig. \. Hot to cold configuration. 


ampoule 


I 


Fig. 2. Cold to hot configuration. 


sample. Fig. 2 depicts the "hot on bottom" config- 
uration, where the SSA vector travels in a direction 
from the cold end to the hot end of the sample. 

To understand which orientation is preferred one 
must first look at the phase diagram [5] of the 
PbSnTe alloy (Fig. 3). The phase diagram summar- 
izes the equilibrium relationship between the tem- 
perature and concentration values of the alloy. 
From Fig. 4a and Fig. 4b. assume that point “a" 
defines the alloy at the molten state (that is, begin- 
ning with a uniform mixture and the sample has 
been melted back to the seed). If one then cools the 
alloy to point "b", in the two phase region, it will 
separate into two parts a solid phase w^hich con- 
tains a lower concentration of SnTe and a liquid 
part which is more concentrated in SnTe. This, in 
fact, is what is meant by "SnTe is rejected into the 
liquid." Further cooling (point "c") results in an 
even higher concentration of SnTe rejected into the 
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Fig. 3. Phase diagram of PbSnTc [5], 



Solkj-Liquicl interface 



Fig. 4. (a) Phase diagram showing the operating lines, (b) Con- 
centration prifilc of PbSnTc in the solid and liquid. 


liquid. In practice, the liquid concentration at the 
solid-liquid interface will increase continuously 
and smoothly along the liquidus line and the cool- 
ing does not take place stepwise. It is shown in 


Fig. 4a as such only for the sake of dramatizing an 
explanation. 

This rich concentration of SnTe just to the right 
of the solid-liquid interface is dramatically different 
than the composition of the specie at the far right 
portion of the liquid. The composition of the hotter 
liquid is predominantly PbTe. PbTe is more dense 
than the SnTe that is near the solid-liquid interface. 
However, the interface is cooler than the bulk 
liquid and it may be noted that the density de- 
creases with temperature [6]. Therefore, this poses 
an interesting question and introduces the primary 
scientific objective for this AADSF experiment. De- 
pending on the orientation of the SSA vector. 
which configuration is more stable for this particu- 
lar class of material? The configuration from Fig. 1 
depicts a solutally unstable (more dense specie^ 
PbTe, on top of sample), thermally stable (cooler 
specie, SnTe, on bottom) configuration. The config- 
uration from Fig. 2 depicts a solutally stable (more 
dense specie. PbTe, on bottom of sample), thermal- 
ly unstable (cooler specie, SnTe, on top) configura- 
tion. Conducting this experiment in the micro- 
gravity environment of space will allow scientists to 
answer this scientific query. 

2.3. Orbital dynamics experiment requirements 

In order to successfully test the scientific ques- 
tions of this experiment it is paramount that the 
SSA vector for both cases be as close as possible to 
being parallel to the growth axis for each sample. 
Recall that the vector for the first case travels in 
a direction from the hot end to the cold region of 
the AADSF (Fig. 1) and the vector for the second 
case travels from the cold end to the hot end 
(Fig. 2). Three primary factors contribute to the 
steady-state acceleration environment at any loca- 
tion away from the spacecraft center of gravity 
(c.g.): drag, gravity-gradient, and rotational effects 
[7]. Drag effects are a function of spacecraft alti- 
tude. attitude, and the time in the 11 -year solar 
cycle that the mission takes place (atmospheric 
density influence). Gravity-gradient effects are 
primarily a function of the distance away from the 
shuttle center of gravity (c.g.) and the shuttle atti- 
tude. Rotational effects are comprised of radial 
contributions (a function of shuttle angular velocity 
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- components, oj. and the distance away from the 

shuttle c.g.) and tangential contributions (a func- 
ght tion of shuttle angular acceleration components, a, 
ent and the distance away from the shuttle c.g). Meticu- 
ght lous orbital dynamics studies using the Simulation 
and Analysis of Multi-Spacecraft On-Orbit (SAM- 
nse SON) 6-Degree of Freedom (6-DOF) model [8], 
ice. developed by engineers from Teledyne Brown En- 
ulk gineering at Marshall Space Flight Center (MSFC), 

de- were conducted prior to the flight to align the SSA 

>ses vectors as parallel as possible to the long axis of the 

ary respective samples (hot to cold, cold to hot config- 

De- urations). Due to the location of the AADSF rela- 
tor, tive to the shuttle c.g., it became evident from the 

cu- studies that it was not possible to provide perfectly 

g- 1 aligned SSA vectors relative to the long axis of the 
cie, samples. A particular nomenclature was used in 
>ler defining how well these vectors were aligned with 
fig* the long axis of the sample. This nomenclature was 
ore based on the ratio of the acceleration component 
aal- along the long axis of the sample, or the “axiaP 
ira- component, to the acceleration component acting 
-TO- perpendicular to the direction of growth, or the 
s to “normal” component. The objective for both ex- 
perimental cases was to obtain as high as possible 
axial/normal ratio conditions during sample pro- 
cessing. 

les- 

the 2.4. Pre-flight and realtime SSA vector results 

c to at the AADSF sample location 

pie. 

s in Orbital dynamics studies of the shuttle were 
1 of completed to determine the attitudes that provided 

ond optimum axial/normal SSA ratios for the AADSF 

end scientific objectives [9]. These shuttle attitudes, 

the along with the predicted axial/normal ratios for the 

>ca- two experiment runs, are provided in Figs. 5-8 

vity respectively. The resulting AADSF attitudes were 

ects forwarded to the flight operations teams at both 

alti- MSFC and the Johnson Space Center (JSC) to 

olar incorporate during realtime AADSF science opera- 

eric tions. 

are The Microgravity Analysis Work Station 
the (MAWS) team, represented by flight control engin- 

uti- eers with extensive orbital dynamics experience, 

dial were located alongside the AADSF science team in 

^city the Payload Operation Control Center (POCC) at 



Fig. 5. AADSF hot to cold attitude configuration. 



Fig. 6. Pre-flight modeled SSA vector at AADSF. hot to cold 
configuration. 



Fig. 7. AADSF cold to hot attitude configuration. 
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Fig. 8. l>re-Hight modeled SSA veclor at AADSF. cold to hot 
configuration. 


MSFC. The MAWS team was responsible for 
monitoring and evaluating the integrity of the SSA 
environment at the AADSF location during sample 
processing. The MAWS team was also responsible 
during the flight for obtaining updated shuttle and 
atmospheric density parameters and. if necessary, 
make attitude modifications based on realtime 
SAMSON 6-DOF simulation runs. Based on real- 


time conditions, the MAWS team recommendc 
a modification to the first attitude to reduce undt 
sirable gravity-gradient effects acting normal to th 
growth axis of the sample. No modifications wer 
made to the second AADSF attitude. Figs. 9 and 1 
provide a snapshot of the realtime MAWS result 
of the SSA environment (as defined in axial/normr 
SSA ratio terms) during realtime sample processin 
of both crystals. From the realtime data in Figs, 
and 10 one can see that it is possible to analytical! 
map the residual acceleration environment on a ve 
hide platform in Low Earth orbit (LEO) in suppoi 
of g-sensitive material science experiments. Real 
time MAWS results from Figs. 9 and 10 show goo< 
correlation with the SAMSON 6-DOF analytica 
solutions from Figs. 6 and 8. From Figs. 9 and H 
one can observe a few characteristics in the data 
From example, observed spikes represent the el 
fects of Vernier Reaction Control System (VRCS 
jet firings that were required to maintain the cor 
rect shuttle orientation during AADSF sample pro 
cessing. Straight slopes without spikes represen 
loss of signal (LOS) periods where the shutth 
is temporarily out of communication with the 
ground. It is interesting to note from Fig. 9 thi 
period around 196 h MET where the axial/norma 
SSA ratio at AADSF decreased significantly. Thi^ 
event was due to the fact that over that time periot 



Mission Elapsed Time (Hours) 

Fig. 9. Realtime snapshot of MAWS results, hot to cold attitude configuration. 
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Fig. 10. Realtime snapshot of MAWS results, cold to hot attitude configuration. 


the attitude control of the shuttle was temporarily 
turned off as a result of thermal problems in the 
VRCS system. The thermal problem was quickly 
resolved and attitude control was reacquired. 
Nevertheless, from this event one observes the sen- 
sitivity of the SSA environment as a function of 
shuttle orientation. 


3. Conclusion 
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ABSTRACT 

The results of a numerical study that models the three dimensional flow of solutal convection in a cylindrical 
geometry is presented. The model was developed to estimate the convective flows that accompany the 
experimental measurement of molecular difflisivity of oxygen ions in liquid metals by the method of 
electrochemical titration. The model predicts convective flows imder certain conditions. The predictions are 
in qualitative agreement with the experimental results that show an enhancement in the effective diffusivity. 

© 1999 COSPAR. Published by Elsevier Science Ltd. 


INTRODUCTION 

This paper is concerned with the numerical modeling of the transport effects during the measurement of 
oxygen diffusivity in liquid metals and the qualitative comparison of the model with experiments. Several 
investigators using electrochemical titrations have measured the diffusivity of oxygen in various liquid 
metals and the procedures are well discussed. (Tare, 1980). The objective of the present investigation is to 
obtain a qualitative picture of the dominant convective processes by employing a three dimensional 
numerical model and also to compare the qualitative results of this model with careful experiments done by. 
Sears et al. (1993) and reconfirmed by the first author of this paper. 

Experiments for diffusivity measurements involve the establishing of an oxygen concentration gradient. 
Calculations that simulate transient and steady state experiments are therefore given for the cases when the 
oxygen is depleted from the liquid metal sample in two different modes- a bottom depletion mode and a top 
depletion one. In each mode the liquid metal sample is assumed to be nearly vertical. In the bottom 
depletion mode the concentration of oxygen at the bottom of the sample is very small and oxygen is depleted 
from that boundary. Therefore the fluid is gravitationally stable in the bottom depletion mode. For this case, 
calculations indicate the existence of stable density gradients whereas in the top depletion mode the reverse 
is true and solutal convection is predicted depending on the experimental conditions. The distribution of the 
oxygen concentration and the flow patterns ^t were obtained from the simulation studies indicate that the 
measured diffusivity can be easily corrupted by solutal convection. This observation is corroborated by 
experimental results. The experimental diffusivity measurements were strongly affected by convective 
effects in those top depletion experiments that had a higher initial oxygen concentration. However, the 
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unsteady state top depletion mode was preferred when the initial oxygen levels were low and the numerica. 
model also confirms this. As evidenced from the modeling, the strength of convection observed in a stead\ 
state top heavy arrangement was much larger than the bottom heavy configuration. The results of botf 
experiments as well as modeling indicate that solutal convection is greatly enhanced depending on the 
experimental conditions. The parameters, which affect the results, are the mode of operation, initial oxygen 
concentration, tilts and heat leakages. 


THE MODEL 

Figure 1 is a schematic of the problem that was modeled numerically to obtain three dimensional flow fields 
The model that was used assumed that the Boussinesq equations hold. The concentration conditions imposed 
on the container depended on the mode and place of oxygen removal. These are shown in Figure 1. The 
calculations were done assuming that the fluid is incompressible and Newtonian. The thermophysical 
properties of tin used in the calculations are well established and are given in the sources cited by Sears et.al. 
(1990). Kao (1993) estimated the solutal expansion coefficient for oxygen in tin as 0.865 /mole fraction, by 
using a dilute solution assumption and a simple mixing rule. The value of diffusivity used in the calculations 
that predict the flow was obtained from Sears, et. al and was assumed to be 5.7 x 10’’ m' /s No-slip 
conditions were used at all boundaries. The governing equations, in the general vector form, for the 
calculations are: 


V-K = 0 


( 1 ) 


5F - - 1 

— 4- F . VF = --VP + vV^F + g[;0,(r - r„) + y5,(c - 


( 2 ) 


- + V-VT 
dt 



( 3 ) 


-.F.VC = DV^C (4) 

where F , P, T and C are the velocity, pressure temperature and concentration fields. The subscript ‘ 0’ 
means a reference state such as the initial state, v , k , p , C ^ and D are the kinematic viscosity, thermal 

conductivity, density, specific heat and mass diffusivity respectively while fJj and /?, are the thermal and 
solutal expansion coefficients. 


The equations were solved numerically using a SIMPLE algorithm, a method developed by Patankar 
(1980). The diameter of the cylindrical cell was assumed to be 0.73 cm as this was the dimension of the 
cell used in the experiments of Sears et. al. and in the present study. The heights of the cell used in the 
calculation varied according to the different cases studied and are reported in Table 1 . 
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Fig. 1 . Schematic representation of the cell used 
in the numerical calculations. 

RESULTS AND DISCUSSION 

As the emphasis in this study was an understanding of the qualitative features of the numerical simulations 
of the physical experiment, these features are discussed along with the experimental results. Figure 2 
represents the unsteady convection that arises in the top and bottom depletion modes during the 
electrochemical titrations.. Note that in both modes of oxygen depletion the container is assumed to be 
slightly tilted (z-axis was 2° from the vertical line) with respect to gravity. This is within reason for it is 
impossible to restrict the container's orientation within 2° in the actual experiments. Also a heat leakage of 2 
degrees is assumed to exist at the sidewalls. At first sight it would appear that the bottom depletion mode 
ought to be solutally stable, as the arrangement is bottom heavy with respect to density gradients. However, 
convection can always be expected on account of the small tilt. Even though the tilt is extremely small, it is 
not negligible causing the concentration gradient to become oblique with respect to the gravity vector and 
generating convection immediately. The top depletion mode on the other hand is conditionally unstable even 
if the container geometry were vertical as the fluid is in a top-heavy arrangement and the tendency to give 
rise to convection is a function of the height of the sample and the initial concentration of oxygen in the tin. 
In fact it may be noted from Figure 2 that the convective flows are strong in the top depletion case. The fact 
that the oxygen is quickly depleted in the top depletion mode is also very apparent and this causes an 
apparent increase in the diffusivity values. This is also seen in the experimental results given in Table 1 . The 
effect of initial concentration is seen in the comparison of Figure 2 with Figure 3. The essential difference 
between these two figures is the initial concentration at the start of the titration procedure. There is however 
a difference in the liquid aspect ratios between Figures 2 and 3. The liquid depths and initial conditions as 
well as mean temperatures that were assumed in order to obtain Figures 2 through 4 correspond to the 
experimental values given in Table 1 while the diameter was assumed to be 0.73 cm. The different liquid 



1306 


C. Miillika et ai. 


aspect ratios are of little consequence as the diffusion boundary layer length is comparable in bot 
calculations. The diffusion boundary layer changes with slowly time and as the effective aspect ratio i 
determined by the diffusion boundary layer, it can be considered constant for all practical purposes. T 
understand how initial concentration plays a role, we recall the unsteady Rayleigh- Benard problem. We the 
observe that under the top depletion scenario a large initial concentration results in a steep unstable densit 
gradient and this gives way to an ever more top-heavy arrangement. Here the large concentration gradier 
only makes matters worse and promotes convection ultimately corrupting the diffiisivity measurements. O 
the other hand when we consider the bottom depletion case we observe that the concentration gradier 
corresponds to a bottom heavy arrangement and that a steeper concentration gradient simply m^es th 
system even more stable. 


Table 1 . Experimental Conditions and Results for the Diffiisivity of Oxygen in Liquid Tin 


Titration mode 

Height of 
tin (cm) 

Temperature 

(K) 

Initial oxygen 
mole fraction 

Diff. Coefficient 
(m^ /s) 

Transient ^ 

Bottom 

0.389 

994 

2.4 X lO"* 

7.57 X 10'^ 

Top 

0.389 

994 

2.4 x 10"* 

2.26 X 10'* 

Transient 

Bottom 

0.520 

845 

7.7x lO ^’ 

3.43 X 10‘® 

Top 

0.520 

845 

7.7 X lO'* 

4.63 X 10'® 

Steady 

state 

Bottom 

heavy 

0.579 

971 

’’6.3 X 10'^ 

3.05 X 10‘* 

top 

heavy 

0.579 

971 

%.3 X 10'^ 

4.98 X 10'* 


Results reported by Sears et. al. (1993) 

The oxygen concentration at the top , The oxygen concentration at the bottom 

From the above it is concluded that large initial concentrations of oxygen in the liquid metal favor a botton 
depletion titration method while for small initial concentrations a top depletion mode can be stable. Figure 
shows that the convection is weak in the top depletion mode because the initial concentration is lower hen 
than in Figure 2. Neither figure shows the concentrations at the depleting surface, which may be assumea 
for all practical purposes, to be zero mole fi'action. The simulations of the steady state transport as shown ii 
Figure 4 depict the convection to be greater here than in the top depletion mode of unsteady transport. Thi. 
happens because the characteristic length in this model is the whole sample height whereas in the transien 
mode it is the thickness of the concentration boundary layer, which in turn is a function of time. Note that thi 
Rayleigh number, a dimensionless group that arises from the scaling of such problems and one which i; 
indicative of the strength of the convection is proportional to the cube of the characteristic length scale (cl 
Narayanan. 1984). It may be noted that the solutal Rayleigh number for the top-heavy system is estimated tf 
be 13859 whereas the calculated value for the bottom heavy configuration is -22595. Hence the strength o 
convection is expected to be more in the top-heavy arrangement which is in accordance with the observation 
However a slight tilt with the steady state mode is far worse than a similar tilt on the unsteady mode because 
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Fig. 2. Comparison of concentration distribution and flow patterns at higher initial oxygen concentration 
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the characteristic length scales are so different and the steady state convection that results because of 
imperfections due to tilts etc. corrupt the diffusivity measurements more severely. 

We conclude that electrochemical titration procedures for diffusivity measurements must be interpreted very 
carefully for the titration mode affects the measured difflisivities in varying degrees. The convective flows 
that occur are influenced by tilts, heat leakages, initial concentration of oxygen and modes of operation be 
they steady or transient. 
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Abstract 

In this study, pattern formation at the bifurcation point from the quiescent state for free 
surface convection in circular containers is determined. A linearized instability calculation that 
employs three-dimensional disturbances in the presence of physically realistic side wall boundary 
conditions is made. The results of the present study will provide a very useful asymptotic limit for 
nonlinear numerical computations. Under restrictive conditions the current calculations check 
favorably with those of earlier workers. The results of these studies are shown to be compatible with 
the experiments of Koschmieder and Prahl. 

Introduction: 

The onset of flow from a quiescent state, when a layer of fluid is heated from below with an 
upper free surface, is one of the classical problems of fluid mechanical instabilities (1) and also the 
subject of many investigations. In this problem, a layer of fluid with an upper free surface is subjected 
to a temperature gradient that is perpendicular to the interface and one of the goals is to find the 
critical conditions and the associated pattern for the onset of flow from an erstwhile quiescent state. 
Curiously, there is a dearth of good experiments on this problem and most of these have been 
conducted in geometries of large lateral extent while the theories that model these experiments 
assume the absence of vertical side walls. When a free surface is absent, convection is driven solely 
by the gravitational field ^ving rise to 'Rayleigh' convection while the inclusion of a free surface 
means that surface tension gradients also determine the onset conditions as well as the planform of 
convective flow. The convective flow driven solely by surface tension gradients is known as the 
Marangoni effect (2) and the physical nature of this problem has been recently reviewed by 
Koschmieder (3). The present paper is concerned with the effect of lateral side walls on convection 
with a free surface. The objectives of this work were to 1) delineate the critical conditions and 
patterns at the onset of free surface convection and indicate, for particular fluid systems, the regions 
where experiments are best conducted and 2) compare the calculations to recent experiments of 
Koschmieder and Prahl (4). 

Instability of the quiescent state is often analyzed by linearizing the boundary value problem 
about the trivial state and inspecting the eigen-spectrum. The associated Marangoni and Rayleigh 
numben are the critical operating parameters or bifincation points while the eigenfunctions represent 
the patterns at the onset point. Theoretical work on convective onset for purely Marangoni flow in 
bounded containers with realistic boundary conditions was done by Vrentas et al. (5), Winters et al. 
(6), Dijkstra (?X Duh (8), Van der Vooren and Dijkstra (9). These papers were, however, restricted 
to calculations of two dimensional disturbances in either a right circular cylinder or in a rectangle. 
Chen et al. (10) attempted to study the pattern formation at the onset in a right circular cylinder but 
the results of their study are suspect in that the large aspect ratio results were not consistent with the 
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wide geometry calculations of the classical theory by Nield (11), at least for the case of 'zero' 
azimuthal dependence. Moreover, the side wall conditions of 'insulation' are rarely compatible with 
experimental conditions that attempt to satisfy the so called 'conducting' case i.e, the situation where 
the side wails have the same thermal conductivity as the operating fluid in question. Recently, 
Wagner, Narayanan and Friedrich (12), assuming a flat surface, have completed three-dimensional 
nonlinear calculations. Here too, the restrictions on the problem were severe in that the gravitational 
field was entirely eliminated and the side walls were assumed to be thermally insulating. Nevertheless 
thdr results were curious in that they indicated a lack of transcritical behavior for small aspect ratios 
and also showed a transition fi’om three to two dimensional flow as the supercritical Marangoni 
number was increased. Note that their results were obtained for the pure Marangoni problem. Other 
studies in bounded geometries for lateral boundary conditions of vanishing vertical component of 
vorticity (cylindrical geometry) or vanishing stress (rectangular geometry) were performed by 
Rosenblat et ai. (13 and 14) and McTaggart (15). Dijkstra (16) recently has given the results for a 
realistic rectangular geometry. A common result of all of the three-dimensional calculations is that 
the aspect ratio and geometry strongly affect the flow structure at and beyond the critical or onset 
point. 

On account of the restrictions either on the dimensionality or the nature of the boundary 
conditions or the gravitational level, these earlier calculations have limitations and are not quite 
compatible with experiments. The present work is motivated by the absence of any 3-D calculations 
of the convection problem in a laterally bounded container where both Raylei^ and Marangoni 
effects are taken into account and also where realistic conditions on the side walls are imposed. We 
expect that such a study with both buoyancy and interfacial tension gradients will provide very useful 
information on the development of planform structure and will provide the real basis for future 
nonlinear numerical studies wWch in turn are necessary for the correct interpretation of experiments. 
This, therefore, is the underling reason for the present communication. 

The Model 

The model that was analyzed describes a liquid such as oil imderiying a gas such as helium 
or air as depicted in Figure 1 . The governing equatioru were derived in a manner similar to that of 
Vrentas et al. (5) with different scale factors. The lower solid boundary has a constant temperature 
Tb and the upper solid surfece has a temperature T, respectively. Tj is the temperature of the interface 
in the conductive state. To simulate the convection the Boussinesq form of the continuity, Navier- 
Stokes and energy equations were used. In what follows the fi'ee surface is assumed to be non 
deformable. This assumption is mildly restrictive for the case of a liquid superposed by a quiescent 
gas as seen by the results of Zhao et al.(17). However in the case of the pure Marangoni problem 
relaxation of this assumption is known to lead to a long wave length instability in an unbounded 
container. The equations of continuity, momentum, and energy were non-dimensionalized using 
different scales (18) and then linear instability theory was applied. The solution technique rests on a 
spectral decomposition in the manner used by Chen et al. (10) and Hardin et al. (19). The 'radial' 
spectral functions that were used depend on the side wall conditions. The flow pattern at the onset 
of convection was obtained fi'om the calculation. These were determined fi'om the eigenfunctions 
while the eigenvalues were identified as the critical Marangoni numbers. The Rayleigh and Biot 
numbers as well as aspect ratios and azimuthal mode 'm* were fixed in any given calculation. Only 
the steady equations were considered because the principle of exchange of stability was assumed. 
Earlier calculations of Vidal and Acrivos (20) as well as Takashima (21) indicate this for the 
unbounded layer case. 
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Discussion of the Results 


General Comments: 

In order to check the reliability of the numerical scheme and the accuracy of the spectral 
representation, the first three critical Marangoni numbers or bifurcation points were calculated for 
the axisymmetric case ( m=0) and Ra=0 for a variety of aspect ratios at Biot numbers equal to 1 and 
100 and these were compared with the corresponding results of Vrentas et al. (5). All of these results 
are depicted in Table 1. Moreover, the first bifurcation point was compared with the corresponding 
results of Vrentas et al. (5) for the case of Ma=0. and Bi=100. and these are shown at the bottom of 
Table 1. We see a discrepancy for the pure Marangoni problem only at small aspect ratios while the 
comparison for the pure Rayleigh problem is very good. This discrepancy is reduced for large aspect 
ratios . 
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Table 1. Comparison of the first three critical Ma for Ra=0 and the first critical Ra for Ma=0 
for various A and Bi assuming insulating side walls and m=0, with Vrentas et aL (5). 


Comparison with Experiments: 

As mentioned earlier, there are very few experiments in bounded geometries for the combined 
Rayleigh-Marangoni problem. Notable are those of Koschmieder and Prahl (4). It is observed that 
in an experiment one cannot fix the Rayleigh number and measure a critical Marangoni number, as 
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the real operating variable is the temperature difference and this variable occurs in both of the 
dimensionless groups. Therefore the critical temperature differences were calculated and the 
corresponding critical Rayleigh and Marangoni numbers were reported in Table 2 for the same 
conditions as reported by Koschmieder and Prahl (4). 


r Case 

m 

Calculated Ma 

Commenta 



0 

109.94 

A-2.16, Ra/Ma-lJ2S. 

Ma-87J, 

n-0 


I 

116,67 

BN 0.85, 

1 

2 

114.83 

depth of air layer “ 0.5 mm 

AT- 7.16 “C 


3 

119.66 

depth of the fluid ■ 1593 mm 


I 

0 

104.50 

A- 2.655, Ra/Ma-0.81, 

Ma-81, 


1 

108.00 

BN 0.696, 

m*2 

1 ^ 

2 

104.70 

depth of air layer ■ 0.5 nun 

AT-8.14«C 


3 

106.70 

depth of the fluid ■ 1109 mm 



0 

98.80 

A- 3.295, IU/Ma-0J26, 

Ma-76J, 


1 

103.10 

BN 0.56, 

m"3 

3 

2 

99J1 

depth of air layer ■ OJ nuB 

AT-9J5*C 


3 

99.81 

depth of the fluid * L699 Bua 



0 

97.90 

A- 4.145, Ra/Ma-0.81 

MaF74, 


1 

100.84 

BN0l7, 

bt-3 

4 

2 

98 JO 

depth of air layer - OJ nun 

AT-7.4«C 


3 

98.32 

depth of the fluid - 2.12 nun 

1 


0 

94.07 

A-2J2, Ra/Ma-3J5, 

Ma-75, 1 


1 

98.80 

BN 1 J9 , 

1 

a 

S 

2 

93.79 

depth of air layer - OJ BUB 

AT-3.77*C 


3 

94.90 

depth of the fluid - 4.22 nua 



Table 2. Comparison of the experimental results of Koschmieder and Prahl with calculations 
for various m and with conducting side wall conditions. 

The thermophysical properties of the silicone oil/air system used in the experiments are given in Table 
3 along with the properties of helium which is a gas that will be used by us in future experiments. The 
results for various 'm' are also given in Table 2. The side wall conditions in the calculations were 
assumed to be of the conducting' type. It is seen that most of our calculated critical Marangoni 
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numbers are above the experimental measurements by about 22%. This therefore requires an 
explanation. 

One reason that the predicted critical Marangoni numbers are uniformly greater than the 
measured Marangoni numbers is partly due to the transcritical nature of Rayleigh Marangoni 
convectioa Recent unpublished nonlinear calculations done by the second author of this study show 
that the transcritical region for the experiments are within 5% of the critical. 

Another reason for the discrepancy is due to possible imperfections induced by side walls. In 
studies involving the influence of side walls, small depths mean that the diameter of the containers 
must also be small. This leads to the possibility that side wall imperfections will influence the onset 
conditions and patterns. These imperfections become more pronounced when the temperature 
differences are substantial and it is to be noted that the critic^ temperature difference for onset 
increases with a decrease in liquid depth. Thus even small imperfections can cause us to lose the 
bifurcation or sudden onset nature of the problem. A weak flow ensues and an increase in the 
temperature difference causes the solution to move along the vicinity of the first branch emanating 
from the first bifurcation point. 

A third reason for the discrepancy between the model and experiments is due to an uncertainty 
in thermophysical properties. The 'over prediction' of the critical Marangoni numbers is entirely 
consistent with the fact that visual detection of the onset point is known within 10% and most but 
not all of the thermophysical properties are known within 10%. The greatest imcertainty in 
thermophysical properties is in the dynamic viscosity. It has been our past e 3 q;)erieoce (22) th»^t Dow 
Coming silicone oils which are labeled 100 centipoise often have a mean dy nami c viscosity of about 
1 S-20% less than what that label would indicate. Hus fiict alone is enou^ to explain the apparent 
discrepancy between our predicted results and those rqKXted by Koscfamieder and Prahl making their 
adjusted results very close to our predictions. These ftcts in ootyunctkm wth the inq>erfectkMis cited 
above can explain the theoretical over prediction of the experinoentaOy determined critical Marangoni 
numbers. 



m 

V 

cmVs 

K 

cmVs 

k 

W/cmK 

WM 

(-do/dT) oil-fluid 
dync/cm. *C 

1 Silicone Oil 

0.00096 

1 

0.001095 

0.001588 

0.968 

0JI5 

1 Air 

0.00333 

0.157 

0.1818 

0.000262 

•4WU 



— 



■iiil 


Table 3. Physical properties of fluids used in typical experiments. 


We now comment on the predicted and observed patterns. The predicted pattern at the onset 
of convection is the same as obtained by Koschmieder and Prahl for the aspect ratio of 2.16 but 
differs from the experimental ones for the other aspect ratios that are reported in Table 2. The 
calculated bifurcation points or critical Ma, that correspond to the various values of *m', are clustered 
for all aspect ratios greater than 2. 16. In fact for the larger aspect ratios it is seen that the first two 
critical Ma are within about 1 % of each other. Now it may be noted firom the work of Tavantzis et- 
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al. (23) and Matkowsky and Reiss (24) that characteristics of the solutions along the first branch are 
strongly influenced not only by the eigenvector at the first bifurcation point but also by the 
eigenvectors that are associated with the subsequent nearby eigenvalues. As a result small 
imperfections on problems with moderate aspect ratios, where bunching of solutions takes place will 
easily cause us to experimentally obtain patterns that are different from the ones that are predicted 
by the first bifurcation point obtained by linear theory which assumes perfect or ideal conditions. In 
other words, mild imperfections and clustering of bifurcation points cause us to believe that there 
is no substantial contradiction between the predicted patterns in our results and those of Koschmieder 
and Prahl (4). While this is so, it suffices to say that experiments are better predicted by calculations 
that consider both the non ideal or imperfect' nature of walls and the transient nonlinear interactions 
so as to allow us to go beyond the onset point of convective flow and also consider transcritical 
behavior. 

Figures 2a through 5b represent the three-dimensional patterns and the corresponding 
planforms at a specific 'z* level for various values of ”m'. It is clear that for 'm' = 0 we have 
axisymmetiy while for 'm* =1 the flow is exactly antisymmetric every z radians. When 'm'=2,3 the 
antisjrometry occurs every n/2 and n/3 radians. The ctilculations are ^own for case 5 in Table 2 and 
it may be seen that the experiments of Koschmieder and Prahl (4) indicate a 'm' =3 flow while 
calculations predict 'm'=2 to be the most unstable. However the modes *m*=0,2, and 3 are within 1% 
of each other and ^ven the difficulty in controlling temperature differences, it is conceivable that 
pattern switching could easily take place. 

We point out that the second author has, in collaboration, with others recently conducted 
experiments in three different aspect ratios ( 2.53, 1.49 and 0.75). The modes that are predicted by 
our method for these aspect ratios are m=0, m=0 and m^l respectively and coincide precisely with 
those observed in these experiments. Moreover the mode mr2 is veiy close to the m=0 mode for the 
aspect ratio 2.53 and the experiments also bear this out by generating pattern switching behavior. 

Summary; 

The critical conditions and patterns for the onset of three-dimensional convection in bounded 
circular containers were obtained. It was found that for a large range of aspect ratios the 
axisymmetric mode or 'nvO* mode is the most unstable and that codimension 2 points occur for either 
large or somewhat small (< 1 .0) aspect ratios. The comparison with existing experimental patterns 
is within reason but not exact and we have rationalized that this is possibly due to the effect of 
inq)er{ectk>ns on geometries that encourage closely bunched dgenmodes. A reason for the apparent 
discrepancy between predicted and experimentally determined critical Ma is that the experimental 
results have assumed the correctness of the viscosity of the test fluid as stated by the manufacturer 
and we believe that the actual viscosity is substantiaUy less. We believe that this is the major reason 
that might explain the difference in reported critical values of the experiment and those predicted by 
the theory. Calculations in deep layers and narrow containers show promise of future experiment^ 
verification. 
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In this study, pattern formation at the bifurcation point from the 
quiescent state for free surface convection in circular containers is 
determined. A linearized instability calculation that employs three- 
dimensional disturbances in the presence of physically realistic 
side wail boundary conditions is made. The results of the present 
study will provide a very useful asymptotic limit for nonlinear 
numerical computations. Under restrictive conditions the current 
calculations check favorably with those of earlier workers. The 
results of these studies are compared with the experiments of 
Koschmieder and Prahl (J. Fluid Mech, 251, 571, 1990). o i996 

Academic Pram, Inc. 

Key Words: Marangoni-B6nard convection; instability; pattern 
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INTRODUCTION 

The onset of flow from a quiescent state, when a layer of 
fluid is heated from below with an upper free surface, is one 
of the classical problems of fluid mechanical instabilities 
( 1 ) and also the subject of many investigations. In this prob- 
lem, a layer of fluid with an upper free surface is subjected 
to a temperature gradient that is perpendicular to the inter- 
face and one of the goals is to find the critical conditions 
and the associated pattern for the onset of flow from an 
erstwhile quiescent state. Curiously, there is a dearth of good 
experiments on this problem and most of these have been 
conducted in geometries of large lateral extent, while the 
theories that model these experiments assume the absence 
of vertical side walls. When a free surface is absent, convec- 
tion is driven solely by the gravitational field giving rise to 
“Rayleigh’’ convection, while the inclusion of a free surface 
means that surface tension gradients also determine the onset 
conditions as well as the planform of convective flow. The 
convective flow driven solely by surface tension gradients 
is known as the Marangoni effect (2) and the physical nature 
of this problem has been recently reviewed by Koschmieder 
( 3 ) . The present paper is concerned with the effect of lateral 
side walls on convection with a free surface. The problem 

' To whom correspondence should be addressed. 


of pattern formation as affected by side walls is strongly 
connected to crystal growth by the Bridgman technique. In 
fact, during crystal growth, the aspect ratio is ever changing 
on account of the moving interface. But this is not the only 
reason why this problem is being studied. Convection experi- 
ments are necessarily done in closed containers and the cor- 
rect interpretation of these experiments requires a theory that 
is closely connected to them. Moreover detailed numerical 
calculations that examine the splitting of solutions and gener- 
ation of cascaded bifurcation as the aspect ratio of the con- 
tainers increase must be checked with appropriate asymp- 
totic models of which the current study is an important one. 

Experiments by Koschmieder and Prahl (4) on convection 
in a laterally bounded container show the marked effect 
that the geometry has on pattern selection. The theoretical 
investigations of the onset of free surface convection in 
bounded containers are restricted and are approximate in 
many ways and none thus far have been posed so as to 
closely replicate experimental conditions. As a result, it has 
not been heretofore possible to examine whether the essential 
qualitative features of the actual problem have been retained 
by the approximate models. It is therefore the objective of 
this work to delineate the critical conditions and patterns at 
the onset of free surface convection and indicate, for particu- 
lar fluid systems, the regions where experiments are best 
conducted. Another objective is to compare the calculations 
to recent experiments of Koschmieder and Prahl (4). 

Instability of the quiescent state is often analyzed by lin- 
earizing the boundary value problem about the trivial state 
and inspecting the eigenspectrum. The associated Marangoni 
and Rayleigh numbers are the critical operating parameters 
or bifurcation points, while the eigenfunctions represent the 
patterns at the onset point. Theoretical work on convective 
onset for purely Marangoni flow in bounded containers with 
realistic boundary conditions was done by Vrentas et ai (5 ), 
Winters et ai (6), Dijkstra (7), Duh (8), and Van der 
Vooren and Dijkstra (9). These papers were, however, re- 
stricted to calculations of two-dimensional disturbances ei- 
ther in a right circular cylinder or in a rectangle. Some of 
these studies indicated the existence of a transcritical region 
in the vicinity of the onset of flow. This means that convec- 
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tive flow could begin ahead of the theoretically predicted 
bifurcation point or onset point. The same studies also 
showed that the transcritical region decreases as the width 
of the container increases. Chen et al (10) attempted to 
study the pattern formation at the onset in a right circular 
cylinder but the results of their study are suspect in that the 
large aspect ratio results were not consistent with the wide 
geometry calculations of the classical theory by Nield ( 1 1 )» 
at least for the case of “zero” azimuthal dependence. More- 
over, the side wall conditions of “insulation” are rarely 
compatible with experimental conditions that attempt to sat- 
isfy the so called “conducting” case, i.e., the situation where 
the side walls have the same thermal conductivity as the 
operating fluid in question. Recently, Wagner et al (12), 
assuming a flat surface, have completed three-dimensional 
nonlinear calculations. Here too, the restrictions on the prob- 
lem were severe in that the gravitational field was entirely 
eliminated and the side walls were assumed to be thermally 
insulating. Nevertheless their results were curious in that 
they indicated a lack of transcritical behavior for small aspect 
ratios and also showed a transition from three- to two-dimen- 
sional flow as the supercritical Marangoni number was in- 
creased. Note that their results were obtained for the pure 
Marangoni problem. Other studies in bounded geometries 
for lateral boundary conditions of vanishing vertical compo- 
nent of vorticity (cylindrical geometry) or vanishing stress 
(rectangular geometry) were performed by Rosenblat et al 
(13, 14) and McTaggart (15). Dijkstra (16) recently has 
given the results for a realistic rectangular geometry. A com- 
mon result of all of the three-dimensional calculations is 
that the aspect ratio and geometry strongly affect the flow 
structure at and beyond the critical or onset point. 

On account of the restrictions either on the dimensionality 
or on the nature of the boundary conditions or the gravita- 
tional level, these earlier calculations have limitations and 
are not quite compatible with experiments. The present work 
is motivated by the absence of any 3-D calculations of the 
convection problem in a laterally bounded container where 
both Rayleigh and Marangoni effects are taken into account 
and .also where realistic conditions on the side walls are 
imposed. We expect that such a study with both buoyancy 
and interfacial tension gradients will provide very useful 
information on the development of planform structure and 
will provide the real basis for future nonlinear numerical 
studies which in turn are necessary for the correct interpreta- 
tion of experiments. This, therefore, is the underling reason 
for the present communication. 

THE MODEL 

The model that was analyzed describes a liquid such as 
oil underlying a gas such as helium or air as depicted in Fig. 
1. The governing equations were derived in a manner similar 
to that of Vrentas et al (5) with different scale factors and 



FIG. 1. Schematic of the physical system. 


so, for brevity, the intermediate details are dispensed with. 
The lower solid boundary has a constant temperature To and 
the upper solid surface has a temperature T, , respectively. 
T{ is the temperature of the interface in the conductive state. 
To simulate the convection the Boussinesq form of the conti- 
nuity, Navier— Stokes, and energy equations were used. 

In order to get nondimensional forms of the above-men- 
tioned equations, the scale factors for radial and axial dis- 
tance, velocity, time, and pressure were introduced. These 
are /?, //, axPHIfx, and cti/?, respectively. Here ctj , 

which is usually positive, is the temperature coefficient of 
surface tension, 0 is the temperature gradient in the liquid, 
K is the thermal diffiisivity of the lower or liquid phase, /x is 
its dynamic viscosity, and H is its depth. The dimensionless 
temperature T was scaled with respect to the temperature 
gradient in the liquid phase. In what follows, the free surface 
is assumed to be nondeformable. This assumption is mildly 
restrictive for the case of a liquid superposed by a quiescent 
gas as seen by the results of Zhao et al (17). However, in 
the case of the pure Marangoni problem relaxation of this 
assumption is known to lead to a long wavelength instability 
in an unbounded container. 

The linearized equations of continuity, momentum, and 
energy will be given below. Here P is the pressure field and 
Ra is the Rayleigh number referred to the lower phase and 
Ma is the Marangoni number, while Pr is the Prandtl number, 
f/, V, and W are the dimensionless radial, azimuthal, and 
vertical components of velocity, respectively. Assuming the 
Boussinesq approximation, the perturbed form for the depen- 
dent variables follows. 

The continuity equation is 


1 5 ^ dV dW 

— ““ ( r f/) + -f- 

Ar dr Ar dd dz 


= 0. 


[ 1 ] 


The r, 6, and z components of the equations of motion 
for an incompressible Newtonian fluid become 
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^ ,dU \ dP 

Pr ^ = TT + 

dt A dr 


- 


_1 
A V 


U - 


1 dV 


aV 


[ 2 ] 


hH 

Bi = — , 
k 


[12] 


Pr 


9V _ 1 dP 


1 




2 dU 


Pr 


91V _ dP 
dt dz 

The energy equation is 


+ Ma*‘Rar. 


[4] 


— = Ma-'V"T -h W. 
dt 


[5] 


Here A is the aspect ratio {RIH) and the Laplacian operator 
is 

^2 ^ d ( d\ 1 9' 9' . 

^ ~Ahdr\dr)^A^r^de^ dz^' 

For the lower solid surface that bounds the system no slip 
was assumed, yielding vanishing velocities. The temperature 
is assumed to be uniform at the bottom. 

At the interface we have Newton's law of cooling 


— -t- Bir = 0. 
dz 


[7] 


A momentum balance at the interface gives the following 
conditions for the tangential component 


dU ^ \ dT ^ dV ^ I dT 
dz A dr dz Ar dO 


[ 8 ] 


where h is the heat transfer coefficient and k is the thermal 
conductivity. 

The thermal conditions on the side walls are either perfect 
insulation or perfect conduction. In the latter case the con- 
ductivity of the side walls are the same as that of the liquid 
so that both will have the same constant vertical temperature 
gradient in the trivial state. 

The trivial solution to the above problem is given by 
the quiescent, conductive state in the liquid layer with the 
hydrostatic pressure gradient balancing the buoyancy. Linear 
instability theory was applied, by expanding the unknown 
variables in a perturbation series around the trivial solution. 
The solution technique rests on a spectral decomposition in 
the manner used by Chen et ai (10) and Hardin et ai (18). 
The “radial” spectral functions that were used depend on 
the side wall conditions. The flow pattern at the onset of 
convection was obtained from the calculation. These were 
determined from the eigenfunctions while the eigenvalues 
were identified as the critical Marangoni numbers. The Ray- 
leigh and Biot numbers as well as aspect ratios and azimuthal 
mode m were fixed in any given calculation. Only the steady 
equations were considered because the principle of exchange 
of stability was assumed. Earlier calculations of Vidal and 
Acrivos (19) as well as Takashima (20) indicate this for 
the unbounded layer case. The three-dimensional linearized 
equations were converted into a simpler set using the trans- 
formations ( 10) 


The equation of state for the interfacial tension was used. 
This takes into account the temperature dependence <7i of 
the interfacial tension in a manner similar to the equation 
of state for the density, in the Boussinesq approximation. 
The Rayleigh and Marangoni numbers are given by 


mr. 

e. 

z) = 

u'ir. 

z)cos 

m0 

V{r, 

0, 

z) = 

V'(r, 

z)sin 

m0 

W(r, 

0, 

z) = 

W’(r 

, z)cos 

; m0 

P(r, 

0, 

z) = 

P'ir, 

z)cos 

m0 

T(r, 

0, 

z) = 

T'(r, 

z)cos 

m0 


[13] 


and 


Ma = 

fJLK 

[9] 

W = 

Ra = , 

[10] 

V = 


d(p 

dz 


fiK 


where a is the coefficient of thermal expansion. 

The Prandtl (Pr) and Biot number (Bi) are given by 


duj 

dz 


1 d m 


[14] 


Pr = 


JL 

PK 


[ 11 ] 


The dependent variables </?, u, and T' were then expanded 
into various spatial components by applying a spectral repre- 
sentation as follows; 
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a 




FIG. 2. (a) Variation of critical Ma number with the number of vertical grids for the case of Ra = 100, Bi - 1, and A - 8. (b) Vanation of critical 
Ma number with the number of radial terms for the case of Ra = 100, Bi = I, and A = 8. 


I Af(z)X;(r) 

y-i 

o; = i BT(z)YTir) 

A 

r' = S Cf(z)Z”(r). [15] 

The Appendix contains some of the essential details on 
the form of the spectral functions. The residuals obtained as 
a result of incorporating the spectral representation were 
made orthogonal to the basis functions as required in a Galer- 
kin technique. This resulted in a set of ordinary differential 
equations that depended on the vertical coordinate. These 
were then finite differenced in the “vertical direction” with 
second-order accuracy and a set of linear homogeneous alge- 
braic equations were obtained where Ra and Bi as well as 
aspect ratio and m were treated as input variables and the 
critical Marangoni numbers were obtained as eigenvalues. 
This last procedure was accomplished using an IMSL code 
G8CRG. 

Figure 2 shows an example of the variation of critical 
Marangoni numbers with the number of vertical grid points 
and number of radial eigenfunctions used. It was typically 
seen that 25 radial terms and 30 vertical grid points sufficed 
to delineate the physical nature of the problem. 

DISCUSSION OF THE RESULTS 
General Comments 

In order to check the reliability of the numerical scheme 
and the accuracy of the spectral representation, the first three 
critical Marangoni numbers or bifurcation points were calcu- 
lated for the axisymmetric case (m = 0) and Ra = 0 for a 
variety of aspect ratios at Biot numbers equal to 1 and 100 


and these were compared with the corresponding results of 
Vrentas et al. (5). All of these results are depicted in Table 
1. Moreover, the first bifurcation point was compared with 
the corresponding results of Vrentas et al. (5) for the case 
of Ma = 0 and Bi = 100 and these are shown at the bottom 
of Table 1. We see a discrepancy for the pure Marangoni 
problem only at small aspect ratios while the comparison 
for the pure Rayleigh problem is very good. A brief explana- 
tion for this is offered. 

The study of Vrentas et al. (5) involved two separate 
cases. These were the pure Marangoni and the pure Rayleigh 
problems. The eigenfunction in the former was the tempera- 
ture field at the free surface and resulted from the consider- 
ation of a matrix of order n where n represents the number of 
radial terms during the spectral expansion. However, when 
Vrentas ei al. (5) considered the pure Rayleigh problem the 
eigenmatrix was of order where now the number of radial 
and vertical terms in the eigenfunction expansion were equal. 
Had they looked at the combined Rayleigh -Marangoni case 
they would have an order n* maoix as well, no matter what 
value of Ra was chosen. Now, the larger the number of 
radial and vertical terms, the more accurate the result. And 
so we can expect the accuracy of the combined Rayleigh— 
Marangoni problem, in the limit of Ra = 0, to be lower than 
the pure Marangoni problem unless a very large number of 
terms are used. In the present study an eigenmatrix of order 
n(k - 1) results where k is the number of finite difference 
intervals in the vertical direction and n is the number of 
radial terms. Therefore the rate of convergence in our prob- 
lem where Ra = 0 is different and lower than that in the 
pure Marangoni case studied by Vrentas et al. (5). This 
is why we have a discrepancy between our results. This 
discrepancy is reduced for large aspect ratios and as expected 
the comparison is quite favorable when the pure Rayleigh 
problem is considered. 

Table 2 shows the calculated lowest critical Marangoni 
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TABLE 1 

Comparison of the First Three Critical Ma for Ra - 0 (Bifurcation Points) and the First Critical Ra for Ma = 0 
for Various A and Bi Assuming Insulating Side Walls and m = 0, with Vrentas et aL (5) 




A = 

1 


A « 2 


A = 4 


A » 8 


Bi 

This 

work 

Reference 

(4) 

This 

work 

Reference 

(4) 

This 

work 

Reference 

(4) 

This 

work 

Reference 

(4) 

Ma for first 

1 

229.5 

206.0 

133.9 

125.0 

119.0 

120.7 

115.3 

117.2 

bifurcation 

point 

100 

4579.0 

4272.5 

3585.0 

3589.4 

3316.4 

3376.1 

3264.0 

3318.8 

Ma for second 

1 

542.8 

531.6 

184.0 

167.1 

131.8 

128.2 

118.4 

119.4 

bifurcation 

point 

100 

7281.0 

7030.5 

4360,0 

4197.5 

3492.7 

3521.4 

3311.0 

3365,7 

Ma for third 

1 

972.3 

1022.9 

292.3 

282.3 

152.9 

149.0 

124.3 

126.1 

bifurcation 

point 

100 

10060.0 

10051.1 

5285.0 

4948.9 

3813.5 

3863.6 

3391.0 

3446.1 

Ra for first 
bifurcation 
point 

100 

1978.6 

1939.0 

1252.9 

1234.9 

II11.3 

1110.8 

1087.5 

1094.5 


numbers for a variety of aspect ratios and various m at differ- 
ent values of Ra and Bi. The last column of Table 2 gives 
the critical Marangoni numbers obtained from Nield (11) 
for the case of infinitely wide layers. It can be seen that the 
calculated results approach Nield’ s results for aspect ratios 
greater than 8 and are within the asymptotic values of Nield 
( 11 ) by 1.2% for all the values of m assumed. It is further 
noted that the results of Chen et aL ( 10) did not approach 
Nield’s values for the case of m = 0 and are therefore in 
disagreement with our calculations. It is felt that this is a 
consequence of inaccuracies in their calculations and in fact 
this is one reason that caused us to repeat the computations. 
We comment on this further in the Appendix. 


Table 3 shows the comparison in critical Marangoni 
numbers between the case of insulated and conducting 
side walls and in nearly all cases it can be seen that the 
case of conducting side wails leads to higher Marangoni 
numbers, thereby indicating greater stability to distur- 
bances. In this regard the problem is similar to the pure 
Rayleigh problem wherein one can show from self-adjoint 
operators that conducting side walls in a container of arbi- 
trary shape leads to greater stability. A similar proof is not 
available for the combined Rayleigh-Marangoni problem. 
Likewise it can be shown from self-adjoint operators that 
the critical Rayleigh numbers in cylindrical containers 
scales with geometry such that 


TABLE 2 

Critical Marangoni Number for Various m, Bi, and Ra for Insulating Walls 


Critical Marangoni numbers for 


Ra 

Bi 

m 

A » 1 

A - 2 

A - 3 

A - 4 

I 

00 

Nield's 

result 

100 

0 

0 

175.4 

82.3 

74.7 

71.1 

68.2 

68.9 



1 

107.8 

87.2 

73.1 

71.6 

68.7 




2 

145.6 

87.3 

72.2 

70.8 

68.1 




3 

235.2 

89.9 

77.9 

69.6 

68.2 


0 

1 

0 

229.5 

133,9 

121.8 

119.0 

115.3 

116.1 



1 

177.7 

133.6 

127.2 

119.5 

115.9 




2 

203.5 

142.5 

119.0 

118.5 

115.5 




3 

297.8 

143.2 

123.6 

118.2 

115.1 


100 

0.2 

0 

184.3 

90.2 

81.6 

78.4 

75.2 

75.9 



1 

119.3 

93.9 

81.4 

78.5 

76.1 




2 

155.0 

96.7 

78.9 

78.4 

75.0 




3 

246.2 

98.2 

84.7 

76.8 

75.2 
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TABLE 3 

Comparison of Critical Ma for Insulating and Conducting Side Walls and Various Ra, Bi, and in 

Critical Marangoni Numbers for 

Wall 


Ra 

Bi 

m 

condition 

/t = 1 

A = 2 

II 

II 

□o 

II 

1 

100 

0 

Insulating 

Conducting 

219.8 

239.5 

120.7 

121.9 

108.5 

109.5 

105.7 

105.6 

101.8 

101.8 

1 

100 

2 

Insulating 

Conducting 

289.8 

369.3 

130.3 

138.7 

110.6 

110.3 

104.7 

105,6 

101.5 

101.8 


g(RaA^) 

a(lM) 


< 0 . 


[16] 


No such analytical result is available for the combined 
Rayleigh -Marangoni problem because of the non-self-ad- 
joint nature of the linearized equations. Nevertheless, calcu- 
lations were performed to inspect the behavior of critical 
Marangoni numbers with aspect ratios and these are shown 
in Fig. 3. Only for the case of m = 1 do local maxima and 
minima occur, whereas for all other values of m a nonmono- 
tonic variation of Marangoni numbers with aspect ratios 
could not be detected. In fact there are regions where the 
Marangoni number changes very little and thereafter for a 
further increase in aspect ratio a drop in the computed critical 
Marangoni numbers is not seen. 


Comparison with Experiments 

As mentioned earlier, there are very few experiments in 
bounded geometries for the combined Rayleigh -Marangoni 
problem. Notable are those of Koschmieder and Prahl (4) . It 
is observed that in an experiment one cannot fix the Rayleigh 
number and measure a critical Marangoni number, as the 
real operating variable is the temperature difference, and 
this variable occurs in both of the dimensionless groups. 



FIG. 3. Critical Marangoni numbers versus aspect ratio for various 
modes for the case of Ra = 100 and Bi = 1. 


Therefore the critical temperature differences were calcu- 
lated and the corresponding critical Rayleigh and Marangoni 
numbers are reported in Table 5 for the same conditions as 
reported by Koschmieder and Prahl (4). The thermophysical 
properties of the silicone oil /air system used in the experi- 
ments are given in Table 4 along with the properties of 
helium which is a gas that will be used by us in future 
experiments, on account of its high thermal conductivity, 
and for which some calculations are reported in this paper. 
The results for various m are also given in Table 5. The side 
wall conditions in the calculations were assumed to be of 
the “conducting” type. It is seen that most of our calculated 
critical Marangoni numbers are above the experimental mea- 
surements by about 22%. This therefore requires an explana- 
tion. 

One reason that the predicted critical Marangoni numbers 
are uniformly greater than the measured Marangoni numbers 
is partly due to the transcritical nature of Rayleigh Maran- 
goni convection. Recent unpublished nonlinear calculations 
done by the second author of this study show that the trans- 
critical region for the experiments are within 5% of the 
critical. 

Another reason for the discrepancy is due to possible im- 
perfections induced by side walls. In studies involving the 
influence of side walls, small depths mean that the diameter 
of the containers must also be small. This leads to the possi- 
bility that side wail imperfections will influence the onset 
conditions and patterns. These imperfections become more 
pronounced when the temperature differences are substan- 
tial, and it is to be noted that the critical temperature differ- 
ence for onset increases with a decrease in liquid depth. 
Thus even small imperfections can cause us to lose the bifur- 
cation or sudden onset nature of the problem. A weak flow 
ensues and an increase in the temperature difference causes 
the solution to move along the vicinity of the first branch 
emanating from the first bifurcation point. 

A third reason for the discrepancy between the model 
and experiments is due to an uncertainty in thermophysical 
properties. The “overprediction” of the critical Marangoni 
numbers is entirely consistent with the fact that visual detec- 
tion of the onset point is known within 10% and most but 
not all of the thermophysical properties are known within 
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TABLE 4 

Physical Properties of Fluids Used in Typical Experiments 



a rc-') 

u 

(cmVs) 

K (cmVs) 

k (W/cm K) 

p (g/cm^) 

{—da/dT) oil-fluid 
(dyn/cm ®C) 

Silicone oil 
Air 

Helium 

0.00096 

0.00333 

0.003326 

1 

0.157 

1.22 

0.001095 

0.1818 

1.770 

0.001588 

0.000262 

0.0015 

0.968 

0.0012 

0.0001627 

0.05 


10%. The greatest uncertainty in thermophysical properties 
is in the dynamic viscosity. It has been our past experience 
(21) that Dow Coming silicone oils which are labeled 100 
cP often have a mean dynamic viscosity of about 15—20% 
less than what that label would indicate. This fact alone is 
enough to explain the apparent discrepancy between our 
predicted results and those reported by Koschmieder and 
Prahl, making their adjusted results very close to our predic- 
tions. These facts in conjunction with the imperfections cited 
above can explain the theoretical overprediction of the exper- 
imentally determined critical Marangoni numbers. 

We now comment on the predicted and observed patterns. 
The predicted pattern at the onset of convection is the same 
as that obtained by Koschmieder and Prahl for the aspect 
ratio of 2.16 but differs from the experimental ones for the 
other aspect ratios that are reported in Table 5. From Table 
5 one can also see that the calculated bifurcation points or 


critical Ma that correspond to the various values of m, are 
clustered for ail aspect ratios greater than 2.16. In fact for 
the larger aspect ratios it is seen that the first two critical 
Ma are within about 1 % of each other. Now it may be noted 
from the work of Tavantzis et aL (22) and Matkowsky and 
Reiss (23) that characteristics of the solutions along the first 
branch are strongly influenced not only by the eigenvector 
at the first bifurcation point but also by the eigenvectors that 
are associated with the subsequent nearby eigenvalues. In 
other words the solution along the branch that emanates from 
the lowest critical Marangoni number has characteristics that 
depend not only on the eigenvector or pattern that is associ- 
ated with the lowest Ma (bifurcation point) but also on the 
eigenvectors that are associated with the subsequent nearby 
higher bifurcation points. As a result, small imperfections 
on problems with moderate aspect ratios, where bunching of 
solutions takes place, will easily cause us to experimentally 


TABLE 5 

Comparison of the Experimental Results of Koschmieder and Prahl with Calculations for Various m and 

with Conducting Side Wall Conditions 


Case 

m 

Calculated 

Ma 

Comments 

Koschmieder and 
Prahl’s results 

1 

0 

109.94 

A = 2.16. Ra/Ma = 1.225, Bi = 0.85. 

Ma = 87.5, m = 0. 


1 

116.67 

depth of air layer = 0.5 mm. depth of 

AT = 7.16°C 


2 

114.83 

the fluid = 2.593 mm 



3 

119.66 



2 

0 

104.50 

A = 2.655. Ra/Ma = 0.81, Bi = 0.696, 

Ma = 81. m = 2, 


1 

108.00 

depth of air layer = 0.5 mm, depth of 

AT = 8.14“C 


2 

104.70 

the fluid = 2.109 mm 



3 

106.70 



3 

0 

98.80 

A = 3.295, Ra/Ma = 0.526, Bi = 0.56, 

Ma = 76.5. m = 3, 


1 

103.10 

depth of air layer = 0.5 mm, depth of 

AT = 9.55°C 


2 

99.31 

the fluid = 1.699 mm 



3 

99.81 



4 

0 

97.90 

A = 4.145, Ra/Ma = 0.81, Bi = 0.7, 

Ma = 74, m = 3. 


1 

100.84 

depth of air layer = 0.5 mm. depth of 

AT = 7.4*’C 


2 

98.30 

the fluid = 2.12 mm 



3 

98.32 



5 

0 

94.07 

A = 2.82, Ra/Ma = 3.25, Bi = 1.39, 

Ma = 75, »i = 3. 


1 

98.80 

depth of air layer = 0.5 mm, depth of 

AT = 3.77“C 


2 

93.79 

the fluid = 4.22 mm 



3 

94.90 
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FIG. 4. Depiction of 3-D profiles and planforms of the velocity eigensolutions at m = 0 with Ra = 305.85» Ma = 94.07, and A = 2.82. 


obtain patterns that are different from those predicted by 
the first bifurcation point obtained by linear theory which 
assumes perfect or ideal conditions. In other words, mild 
imperfections and clustering of bifurcation points cause us 
to believe that there is no substantial contradiction between 
the predicted patterns in our results and those of Kosch- 
mieder and Prahl (4). While this is so, it suffices to say that 
experiments are better predicted by calculations that consider 
both the nonideal or the “imperfect” nature of walls and 
the transient nonlinear interactions so as to allow us to go 
beyond the onset point of convective flow and also consider 
transcritical behavior. The utility of the present study is to 
provide a good starting point toward such nonlinear calcula- 
tions. 

Figures 4 through 7 represent the three-dimensional pat- 
terns and the corresponding planforms at a specific z level 


for various values oi m. It is clear that for m = 0 we have 
axisymmetry, while for m = 1 the flow is exactly antisym- 
metric every tt radians. When m = 2 or 3 the antisymmetry 
occurs every tt/ 2 and tt/ 3 radians. The calculations are 
shown for case 5 in Table 5 and it may be seen that the 
experiments of Koschmieder and Prahl (4) indicate a m = 
3 flow, while calculations predict m = 2 to be the most 
unstable. However, the modes m = 0, 2, and 3 are within 
1% of each other and given the difficulty in controlling 
temperature differences, it is conceivable that pattern switch- 
ing could easily take place. 

We point out that the second author has, in collaboration, 
with others recently conducted experiments in three different 
aspect ratios (2.53, 1.49, and 0.75). The modes that are 
predicted by our method for these aspect ratios are m = 0, 
m = 0, and m = 1, respectively, and coincide precisely with 


b 




no. 5. Depiction of 3-D proiiies and planforms of the velocity eigensolutions at m = 1 with Ra = 319.35, Ma = 98,28, and A = 2.82. 
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HG. 6. Depiction of 3-D profiles and planforms of the velocity eigensolutions at m = 2 with Ra - 304.85, Ma - 93.8, and A 2.82. 


those observed in these experiments. Moreover the mode m 
= 2 is very close to the m = 0 mode for the aspect ratio 
2.53, and the experiments also bear this out by generating 
pattern switching behavior. 

Other Calculations 

Finally, a set of results are presented in Table 6 for various 
fluid depths and Biot numbers, some of which may be veri- 
fied experimentally. Using the thermophysical properties of 
silicone oil/air and various depths of liquid the critical Mar- 
angoni numbers are given for various values of m. Another 
parameter that is fixed is the Biot number and two values, 
viz. 1/7 and 1, were chosen. The value of Biot number of 
1 n corresponds to an air thickness equal to the liquid depth. 


while the Biot number of unity corresponds to the use of 
helium as the upper gas of the same thickness as the liquid. 
The relative importance of Rayleigh to Marangoni effects 
arc given by the ratio (Ra/Ma = pagH^loy). This group 
is independent of the temperature difference and may be 
adjusted mainly by changing the depth of liquid or the gravi- 
tational level. It represents the importance of buoyancy or 
gravitational effects compared to surface tension gradient 
effects. 

It is observed from Table 6 that, regardless of the value 
of Ra/Ma, the critical Ma are closely bunched for aspect 
ratios greater than 1.5. If one wishes to study the effect of 
side walls and also simulate a low-gravity environment by 
performing ground-based experiments, it is necessary to 
have very small depths and consequently very small diame- 




RG. 7. Depiction of 3-D profiles and planforms of the velocity eigensolutions at m = 3 with Ra - 308.43, Ma - 94.9, and A 2.82. 
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TABLE 6 


Critical Marangoni Numbers for Various Ra/Ma, m, and Bi 


Biot 

Depth 

Ra/Ma 

m 



Critical Marangoni numbers for 



4 = 0.5 

A = 1.0 

4 = 1.5 

4 = 2.0 

4 = 2.5 

4 = 3.0 

1/7 

1 

0.1823 

0 



206.75 

121.49 

100.66 

94.11 

90.06 




1 

— 

214.47 

133.57 

105.58 

95.95 

92.95 




2 

— 

266.29 

147.15 

110.25 

95.32 

89.45 




3 

— 

330.53 

168.64 

118.15 

98.70 

90.84 


2 

0.73 

0 

— 

197.58 

114.82 

95.25 

88.78 

84.88 




1 

— 

205.20 

126.81 

99.80 

90.54 

87.72 




2 

— 

255.89 

139.89 

104.29 

89.96 

84.34 




3 

— 

319,90 

159.48 

111.82 

93.13 

85.67 

1 

1 

0.1823 

0 

— 

247.10 

153.38 

132.12 

125.61 

120.50 




1 

— 

258,30 

169.90 

138.87 

128.90 

125.47 




2 

— 

308,66 

180.84 

141.44 

125.85 

1 19.79 




3 

— 

376.46 

202.74 

148.65 

128.62 

121.30 


2 

0.73 

0 

— 

234.85 

143.80 

123.82 

117.30 

112.45 




1 

— 

245.67 

160.00 

129.97 

120.55 

117.30 




2 

— 

295.20 

170.47 

132.55 

117.65 

112.30 




3 

— 

363.10 

191.70 

139.44 

1120.33 

113.22 

1/7 

5 

4.56 

0 

661.12 

147.30 

81.65 

67.50 

62.97 

60.11 




1 

553.41 

153.90 

91.97 

71.22 

64.38 

62.38 




2 

869.40 

195.30 

101.57 

74.62 

63.86 

59.75 




3 

1104.44 

253.67 

1 17.87 

80.10 

66.27 

60,70 


ters. Naturally this raises the issue of imperfections and 
moreover the small cell diameters make physical visualiza- 
tion of the flow difficult. On the other hand larger aspect 
ratios cause the clustering of eigenmodes and patterns and 
all of the issues raised earlier become relevant. It may also 
be seen from Table 6 that the general nature of our comments 
remain unchanged merely by changing the Biot number. 
Further it is observed from Table 6 that at aspect ratio 2.5 
the most critical mode corresponds to m - 0 and at aspect 
ratio of 3 the most critical mode is m = 2. This means that 
a codimension 2 point occurs between aspect ratios 2.5 and 
3. There is little doubt that such points are difficult to capture 
experimentally, once again because of the clustering of the 
bifurcation points or critical Ma. However, it is also observed 
from the calculation for a depth of 5 mm that another codi- 
mension 2 point (showing the switch between an m = 1 
mode and an m = 0 mode) occurs between aspect ratios 0.5 
and 1. Given these low aspect ratios it is conceivable that 
this multiple point can be captured in a careful experiment 
with reasonably deep layers. 

From our calculations the following general statements 
may be made. First, narrow aspect ratio containers give rise 
to definite patterns at the onset of convection that are stable 
because the bifurcation points are not closely bunched. We 
expect that very large aspect ratio containers will simulate 
the infinitely wide case of classical theory, giving rise to 
hexagons as seen in earlier experiments. However the me- 


dium aspect ratio containers (aspect ratios >3) will show 
spectral crowding and close branches. As it will be difficult 
to control the temperature differences accurately, it will be 
hard to get good agreement between theory and experiment 
at these intermediate aspect ratios. In other words it is better 
to conduct experiments either in very large geometries where 
end effects are negligible or in very narrow aspect ratio 
containers where spectral crowding is absent. While the me- 
dium aspect ratio containers could give theoretically pre- 
dicted Marangoni numbers, the predicted patterns will be 
difficult to obtain experimentally because of the spectral 
crowding. 

Second, experimental studies in narrow aspect ratio con- 
tainers are best conducted in deep layers so that the radius 
of the containers may be large and side wall imperfections 
such as a mismatch in thermal conductivities and a “min- 
iscus” may be minimized. Minscus problems are likely to 
be encountered in narrow containers. In deep layers “ 
Ra/Ma will be necessarily large and if we are to get predicted 
patterns in narrow aspect ratio containers, it will not be 
possible to bias the problem in favor of the Marangoni effect 

^ It is of course possible to choose a liquid with a very low density and 
thermal expansion coefhcient: however, most experiments are conducted 
with silicone oils because of their nearly constant properties and clean 
interfaces. As a result Ra/Ma is typically targe when the depths are greater 
chan 2 mm. 
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and minimize the Rayleigh phenomenon. The effect of all 
of this is that in narrow aspect ratio containers convective 
instability experiments will invariably be buoyancy driven 
with a Marangoni perturbation and not the other way. 

As a final note to readers of this journal we could take a 
look at this paper from a different perspective and use the 
calculations and companion experiments to predict the sur- 
face tension gradient idddT)^ provided that we have a 
good idea of the other thermophysical properties. In other 
words these bifurcation calculations are a good approxima- 
tion to carefully controlled experiments and can be used for 
parameter identification. 

SUMMARY 


u = 

< 

II 

II 

= 0. 

[A4] 

functions 

for mode m = 

= 0 are 


XUr) 

Ma°r) 


[A5] 


/iK) 

/.«) 

YJir) 

= 0 


[A6] 

Z°{r) 

= U6°r), 


[A7] 


where / and J are Bessel functions of the first kind and a° 
and 6° are the roots of 


The critical conditions and patterns for the onset of three- 
dimensional convection in bounded circular containers were 
obtained. It was found that for a large range of aspect ratios 
the axisymmetric mode or m = 0 mode is the most unstable 
and that codimension 2 points occur for either large or some- 
what small (<1.0) aspect ratios. The comparison with ex- 
isting experimental patterns is within reason but not exact 
and we have rationalized that this is possibly due to the 
effect of imperfections on geometries that encourage closely 
bunched eigenmodes. A reason for the apparent discrepancy 
between predicted and experimentally determined critical 
Ma is that the experimental results have assumed the correct- 
ness of the viscosity of the test fluid as stated by the manufac- 
turer and we believe that the actual viscosity is substantially 
less. We believe that this is the major reason that might 
explain the difference in reported critical values of the exper- 
iment and those predicted by the theory. Calculations in 
deep layers and narrow containers show promise of future 
experimental verification. 


+ Ma])h{aJ) = 0 

[A8] 

7,(6;) = 0 

[A9] 

and 


7j(a;)/,(a;) + Ma°)hia]) = 0 

[AlO] 

7o(6“) = 0 

[All] 

for perfectly insulating and perfectly conducting side walls, 
respectively. For m > 1 the trail functions are 

y-., J^iccTr) Ua-r) 

^ ^ 7„(a7) /„«) 

[A12] 

F7(r) = J„iPTr) 

[A13] 

Zfir) = 7„(67r), 

[A14] 

where aj, P”, and 6” are the roots of 



APPENDIX 

The solutions of Eqs. [1] - [5] were sought for the bound- 
ary conditions 

atz = 0 t/=V=W=r = 0 [Al] 

dz A dr dz Ar dd dz 

[A2] 


(«")/„(«;) + y„(a 7 )/„+,«) = o [ai5] 

= 0 [A16] 
- ^..(^7) = 0 [A17] 

/„^,(a7)/„(a7) + 7„(a7)/„+,(a7) = 0 [A18] 

JmidT) = 0 [A19] 
J„{8T) = 0 [A20] 


for perfectly insulting walls (r = 1.0) 

U=V=W = -^=0 [A3] 

A dr 


for perfectly insulating and perfectly conducting side walls, 
respectively. 

As noted in the main text of this paper Eqs. [13] -[15] 
were substituted into the governing equations. The pressure 
was then eliminated and equations of fourth order resulted. 
On using the Galerkin method many integrals result from 
making the residuals equal to zero. An example is 


and for perfectly conducting walls 
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[integral]*; = J r\^D* + ^ D? - — D). 
. (3 + nr) ^ (3 - 3m^) 


- 


XjXtdr. [A21] 


When m = 0, [A5] to [A7] were used and when m = 2, 
3, etc., [A12] to [A17] were used. However, when m - 1 
this and several other integrals gave rise to singularities at 
r = 0. There are at least two ways to get over this problem. 
One way is to look at the geometry of the problem as a 
limiting case of an annular compartment with a very small 
inner cylinder and thereby exclude r = 0 from the domain. 
Another way is to choose the m = 2 spectral functions in 
order to express the flow pattern at m = 1 . This is legitimate 
as the eigenfunctions are complete. We chose to do the latter 
and note that Chen et al ( 10) indicated that [ A12] to [ A14] 
with m - 1 were used in their calculations. As this leads to 
singular behavior, this may be one reason why there is a 
difference between our results and Chen et al. (10). 
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ABSTRACT; 

A 3D numerical calculation is performed on a model that depicts the flow profiles due to thermo-solutal 
convection in a cylindrical tube. The calculations were done with the purpose of delineating the qualitative 
features of the flow profiles for the cases when the container's axis is perfectly aligned with respect to the 
mean gravity vector and also when it changes periodically with respect to the gravity vector. It is found that 
the flow profiles are similar to those of the Rayleigh-Benard problem in the case of perfect alignment while a 
swirling pattern appears when the tube’s axis is not aligned with the gravity vector. This indicates that it 
might be preferable to have a slight tilt in the container axis during crystal growth as swirling flow will 
diminish axial mixing. The solutal convection is the dominant feature of the flow and is affected 
considerablv by the gravity level. © 1999 COSPAR. Published by Elsevier Science Ltd. 


INTRODUCTION 

This is a brief report describing the flow profiles that are induced in a low gravity environment in a 
Bridgman tube in which the fluid occupies a constant volume. The Bridgman tube as considered in this 
study is merely a circular cylinder that is subjected to radial thermal gradients and axial solutal gradients. 
Typically, the Bridgman tube is used in the vertical directional solidification of compound semi conductors 
such as Lead Tin Telluride. The growth of such materials is affected substantially by the convective flow 
profiles that accompany the process. This convection is due to thermal and solutal gradients that are 
generated because of the solidification process. Arnold et. a/. (1991) did calculations to model a GaAs space 
experiment and concluded that three-dimensional flows occur under certain gravitational values and 
orientations. Their calculations were not concerned with solutal convection. Naumann and Baugher ( 1 992) 
have made analytical estimates of radial segregation in Bridgman growth for low-level steady and periodic 
accelerations. In any actual growth process, the liquid zone is ever shrinking and this can be expected to 
change the flow profiles quite a bit. Nonetheless it would be interesting to have an idea of the flow profiles 
that are generated when the force conditions on the ampoule are compatible to a time dependent 
microgravity level and where certain assumptions such as a constant liquid zone is assumed. We present here 
a numerical model that shows the effects of off axis and a time dependent orientation on the flow profiles in 
a Bridgman tube. The effect of tilting the otherwise vertical container with respect to gravity is also 
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described. The calculations were done with the purpose of delineating the qualitative features of the flow 
profiles for the cases when the container's axis is perfectly aligned with respect to the mean gravity vector 
and also when its axis periodically changes with respect to the gravity vector. The gravitational levels that 
are assumed range from 10‘^ge (or ten micro g) to 10 "^gc where ge is earth's gravity. A value often micro g is 
reasonable as it is a fair representation of the low frequency accelerations experienced on the cargo bay of the 
U.S. space shuttle or on the future international space station if the Bridgman tube. A level of 100 micro g is 
not very probable, however we also present calculations that include this extreme case. To the best of our 
knowledge this is the only study that shows the dominant effect of solutal convection over thermal 
convection at the higher gravity levels and also the only study where the ampoule axis orientation is varied 
with time. 


Table 1 The Thermophysical Properties Used in the Calculations 


Density 

7.04 g/cm^ 

Kinematic Viscosity 

0.0024 cm Is 

Thermal Diffusivity 

0.03 cm'/s 

Solutal Diffiisivity 

7 10’^ cm^/s 

Thermal Expansion Coefficient 

1.18 10“* /C° 

Solutal Expansion Coefficient 

0.22 AVeight fraction 

Segregation Coefficient 

0.7 


THE MODEL AND THE NUMERICAL SCHEME 

The model that is used assumes that the Boussinesq equations hold. Further the calculations were done 
assuming that the fluid is Lead Tin Telluride reflecting our interest in compound semiconductors. The 
thermophysical properties of Lead Tin Telluride as used in the calculation are given in Table 1. Figure 1 
describes the situation when a container is subject to thermal gradients with a solidifying interface at 
z=H-Hsoiid. The thermal and concentration boundary conditions imposed on the container are given in Figure 
1. No-slip conditions are used at all boundaries including the solid-liquid interface upon which the 
coordinate system is fixed and which is assumed to move down in the z direction at a constant speed Vs 
equal to 1 cm/hr. The height of the tube of diameter equal to 1. cm.is given by H and assumed to be 5.0 cm., 
equally divided between the solid and liquid zones while the insulation zone is assumed to occupy the 
middle one third. These correspond roughly to the experimental ampoule used in a Lead Tin Telluride 
experiment that was conducted on USMP 3. The hottest temperature is assumed to be 1150 degrees Celsius 
while the coldest temperature is fixed to be 700 degrees, the interface being at 900 degrees. The orientation 
of the gravitational acceleration is expressed in terms of the angle between the gravity and the negative 
direction of the z-axis. The major assumption is that the liquid length is kept constant. As a result it is 
assumed that the end of the liquid was at a constant concentration , Co equal to 0.2 weight fraction. This is 
tantamount to a continuous feeding of such liquid at the solidification rate Vs . Before we go on it might be 
useful to point out that the thermal Rayleigh number is estimated to be about 65 for a gravity level 1 0*^gc 
while the solutal Rayleigh using the same length scale is about 147000 
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Fig. 1 The schematic of the geometry and thermal and 
concentration boundary conditions ofthe calculation. 


RESULTS OF THE CALCULATIONS AND CONCLUSIONS 

The finite volume method. SIMPLE (cf. Patankar,I980) was chosen to solve the governing equations. 
Figure 2 shows the flow profiles at two different gravity levels. Figure 2a describes the pattern that is 
exp>ected at a constant g level of 1 0 ^ge . What is to be observed is vertical stacking of an axisymmetric or 
torroidal pattern. This vertical stacking may be expected as the top of the ampoule is hotter than the bottom 
and the lower cell' is in the insulation zone. The configuration acts like a fluid that is heated from above' 
and the weak flow is primarily driven by radial gradients. The weakness in the lower cell is primarily due to 
the effect of the presence of the no slip' solid boundary. 

The situation changes somewhat for the case of a g level of lO^ge as seen in Figure 2b, for here the solutal 
convection begins to play a part. The solutal gradients are unstable in the sense that they promote flow even 
if the thermal expansion coefficient is negligible. As observed earlier the solutal Rayleigh number is about 
147000 whereas the critical solutal Rayleigh number, in the absence of thermal gradients for this aspect ratio 




3D Numerical Model tor Flow m a Bridgman Tube 


1401 




Fig.3. The 2D flow field near the interface under 10 * . The orientation of the gravity vector is a 

function of time ( a = 1 0“ sin cot , co = 2ti / 3600 ). 
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(radius/ height =0.2) turns out to be about 50000 from the calculations of Hardin et. al.. The unicellular 
patterns that are observed are a result of the soluially driven convection and are expected for this geometry 
according to the results of the Rayleigh-Benard problem (cf. Hardin et.al.. 1990). Notice fiirther that the 
vertical stacking arrangement disappeared for the larger ’g' level indicating the dominance of the solutal 
convection over the thermal convection. Figure 2c shows the flow profiles resulting from a time invariant 
change of the ampoule axis with respect to the mean gravity vector of lO'^ge. What is immediately apparent 
is the swirling flow that helps contain the rejected solute near the solid liquid interface. The 'z' component of 
this type of flow is much smaller than the other two components except near the interface and the end of the 
liquid where all components are set to zero on account of no-slip. It is a concentration induced flow because 
the velocity components near the interface in Figure 2c are much larger than in Figure 2a. In a real crystal 
growth configuration such an off axis tilt would help prevent axial mixing and therefore be beneficial to the 
crystal. If the g’- value is increased by an order of magnitude the flow is mostly of a unicellular style except 
near the interface. This is shown in Figure 2d. Figure 3 depicts the periodic change of the cylinder axis with 
the gravity vector. The frequency was set to be one cycle per hour. This was an arbitrary choice even though 
aerodynamic drag causes a readjustment every 20 minutes or so in a typical space orbiter. All the same the 
results can be expected to be qualitatively similar to those reported here in the case when the frequency is 
increased three fold. The flow profiles at a constant 'z' plane near the solid liquid interface are given at every 
quarter cycle, once a periodic steady state is reached. Notice that the direction of the swirl changes every 
half cycle i.e., when the gravity vector crosses the cylinder axis leading to local mixing near the solid liquid 
interface. 

The effect of gravity level and time periodic off axis alignment show collectively that the convection at low 
gravity in a bottom or top seeded Bridgman tube is primarily in the solutal driven mode as long as the gravity 
level is not very small. Moreover a slight tilt with respect to gravity causes the fluid flow to go into a 
swirling mode so that solute is contained near the solute-generating boundary. 
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The fluid physics of buoyancy-driven (Rayleigh) and interfacial tension-driven 
(Marangoni) convection is examined for two superimposed layers of fluids. This con- 
vection occurs on account of temperature gradients that are imposed perpendicular 
to the fluid fluid interface. Interfacial deflections, small as they may be, play an 
important part in identifying the mechanism that governs the flow, and calculations 
have made that indicate whether hot or cold fluid flows towards or away from a 
crest or a trough. As a result, four possible flow structures or ‘modes' at the interface 
have been identified. Two heating styles, heating from below and above, are com- 
pared and the behaviour of the fluid physics as a function of total fluid depths, depth 
ratios and gravity levels is explained. Changes in modes result because of changes in 
these parameters. We have given plausible physically based arguments that predict 
the sequential change in modes as these parameters are changed and have 'veri- 
fied' our conjectures with calculations. Flow mechanisms in the case of a solidifying 
lower pliase have also been studied, as this has an application to liquid-encapsulated 
crystal growth. Where convection is deemed detrimental to crystal homogeneity, we 
conclude that the liquid-encapsulated method of crystal growth is best conducted 
under Earth's gravity. 


1 . Introduction 

This paper is concerned with the study of convection in fluid bilayers. Interfacial- 
driven convection must necessarily involve at least two fluid layers and we could w^ell 
imagine that the fluid physics of motion, driven by interfacial tension and density 
gradients, depends largely on the heating direction, fluid depths as well as property 
ratios. One motivation for this study stems from an interest in liquid-encapsulated 
crystal growth where a vertical cylinder with thermally insulated side walls encloses 
the melt. The crystal solid phase can be below the melt phase and this corresponds 
to the bottom seeding situation where the liquid melt is now heated from above. An 
encapsulant is often placed above the melt in order to provide a diffusion barrier to 
high volatile constituents in the melt. An example is the growth of gallium arsenide, 
wherein arsenic is the highly volatile component and boron oxide the encapsulant. 
A bilayer with a common interface is thereby created. As a temperature gradient is 
applied across the interface, there are basically two mechanisms which can generate 
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Figure 1. Schematic of the bilayer system. Dashed line is the Hat interface in the quiescent 

state. 

convection, i.e. buoyancy and inter facial tension gradients. These two mechanisms 
are called Rayleigh and Marangoni effects, respectively. In a model problem we could 
apply a temperature gradient that is either parallel or antiparallel to the gravitational 
field and the configuration represents an instability problem which is associated with 
a bifurcation from the quiescent state to the convective state. By applying a linear 
stability theory, we get the sufficient conditions for the onset of convection as well 
as the most dangerous wavelength of an imposed infinitesimal disturbance. 

We can understand the physical mechanisms which are involved in interfacial ten- 
sion gradient convection by considering a bilayer configuration, as shown in figure 1. 
Let 7"i > Till ^nd further assume that gravitational effects are negligible. Now sup- 
pose we give a perturbation to an erstwhile flat interface so that the temperature at 
the point *c* is higher than at *b'. As most fluid bilayers have a negative interfacial 
tension gradient, the inter facial tension at 'c' will be lower than at ‘b' and fluid is 
driven from ’c^ to 'h\ Fluid from both phases must then rush towards *c' and the 
final steady state will depend on tiie fluid property ratios and heights. If we have 
a liquid gas system where the upper gas phase is assumed to be p^issive, then only 
liquid from below will move towards *c' and it follows tliat unless the temperature 
gradient is reversed, the perturbations must decay. Gravity stabilizes or not accord- 
ing to the heating arrangement. It may be pointed out that the mechanism for flow 
can take place even if the interface is always restricted to be flat b('cause temperature 
perturbations are still allowed. However, our interest will mainly focus on the general 
case where interfacial deformations are included. 


2. Earlier work on bilayer convection 

The first theoretical work in Marangoni convection was by Pearson (1958). wherein 
the liquid layer was assumed to be superimposed by a quiescent gas. One of the 
early studies in bilayers was by Sternling & Scriven (1959), who considered the pure 
Marangoni problem using mass transfer as an analogue to heat transfer, and Smith 
(1966). who examined the case of thermocapillary and gravity waves. Sternling & 
Scriven found that convection or instability can occur if the transfer takes place in 
either direction. This problem was later extended to include surface deflection and 
surface viscosity, but it was then assumed that the upper phase was passive. 

This work was followed by Zeren Reynolds (1972), who determined the critical 
temperature gradient for the onset of convection in a bilayer of water and benzene in 
order to compare theory with experiment. Now, in dimensionless form, the critical 
temperature gradient is represented by either the Marangoni or the Rayleigh number. 
As these groups are related to each other by a factor that contains physical properties, 
it is sufficient to calculate either the critical Rayleigh or critical Marangoni number 
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for the onset of convection. When the bilayer was ‘heated from below', Zeren k 
Reynolds (1972) found that the onset of convection was either a buoyancy-driven 
flow generated from the upper phase or an interfacial tension gradient-driven flow 
that started at the interface. When they considered the Cf\se of a low liquid layer 
depth of the lower layer compared with the thickness of the upper layer, they found 
that tlie Marangoni convection in the bilayer served to delay the onset of motion and 
that the flow was primarily a buoyancy mechanism which was driven from the upper 
layer. This resulted in vertical stacking in the upper layer where the upper cell in 
the upper layer was associated with the buoyancy mechanism. On the other hand, 
when higli depth ratios of the lower phase were considered, the onset of flow was 
due primarily to a Marangoni mechanism and not influenced much by the Rayleigh 
effect . 

Following the work of Zeren k Reynolds (1972). Perm & \Wllkind (1982) per- 
formed detailed calculations for the silicone oiPair system with the hope of compar- 
ing their results with the experiments of Koschmieder (1967) and Palmer & Berg 
(1971). All calculations were performed for the case of the bilayer being ‘heated 
from below'. In an effort to trace where the Marangoni regime was distinct from the 
Rayleigh regime, a series of calculations were performed and plotted as the critical 
temperature gradient against the depth of lower layer. They claimed that a drastic 
change in slope of this curve indicated the depth ratio of the lower layer when one 
mechanism took over the other. 

Recent work on bilayer convection includes the interesting studies of Rasenat et 
ai (1989) and Wahal & Bose (1988). Rasenat et ai (1989) investigated the case of 
negligible surface deflections and uncovered oscillatory behaviour. They also con- 
sidered a case of finite interfacial deflections but with negligible interfacial tension. 
While it may be argued that interfacial deflections are very small in comparison to 
the fluid depths, it is our view that the interfacial morphology helps to identify the 
controlling factors of competing convective mechanisms. This is what we aim to show 
in the subsequent sections. 

It is noted that Zeren & Reynolds (1972) calculated the energy contributions 
from the buoyancy and surface mechanisms, as well as the critical Rayleigh and 
corresponding Marangoni numbers, in an effort to trace the leading characteristics 
of the flow. It is obvious that the energetics of the flow are calculated across the 
entire domain of the flow field and give some useful global information. However, 
it is also possible to consider the local behaviour of the flow at the interface using 
linear stability methods and by evaluating the eigenfunctions. We feel that this leads 
to vital information on the flow mechanism at the interface. In particular, we can. 
as we shall see. decide whether we have hot or cold fluid rushing towards or from a 
crest or a trough. We will observe that w(' can have four possible flow modes at the 
interface. This paper concerns operating parameters and, as we change the gravity 
level, the total depth or depth ratio, a sequential change of flow structures is obtained. 
The generality of this sequence depends on whether we Inive a liquid- g^is system, a 
liquid liquid system and whether the bilayer is ‘heated from below or above*. The 
order of flow structures as we change these parameters also depends on whether the 
temperature coefficients of interfacial tension and density are negative or positive. 
From this we can decide whether the flow mode at the interface is promoted by 
Rayleigh or Marangoni effects. 
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3. Theoretical development 

The model that we analyse is schematically shown in figure 1 . The governing 
equations are derived in a manner identical to that of Perm fc Wollkind (1982) and 
so we dispense with the intermediate details in the cause of brevity. Without loss 
of generality, we introduce a two-dimensional coordinate system. In this system, the 
direction of the gravitational acceleration corresponds to the negative z-direction. 
The position of the deflecting interface is a function of x and the time t and is 
measured from the datum r = 0. The upper fluid is designated with + and the 
lower fluid with — , so that represents the vertical depth of the upper fluid and 
d the depth of the lower fluid. The lower solid -liquid boundary has a constant 
temperature and the upper solid surface has a temperature of T\. respectively. 
The static interface has a temperature of Tj. We will use the Boussinesq form of the 
continuity, Navier-Stokes and energy ecpiations. 

In order to get non-dimensionalized forms of the above equations, we introduce the 
scale factors for distance, velocity, time and pressure. These are K/d“, {d ~)'^ /k 
and jiK/[d~Y . respectively, k and /i are the thermal difFusivity and dynamic viscosity 
of the lower phase, respectively. 

The dimensionless temperature 0 is defined as 


iT-T.) 

(T,, - TO' 


(3.1) 


In what follows, several important dimensionless groups will arise. These are the 
Rayleigh number i?, Marangoni number 4/, Cripsation number C, the Bond number 
G and Fraud tl number P. They are defined as follows: 


^ ^ ^ cry3(P ^ ^ ^ p ^ 

KV Kfl aod <T() K 

All the physical properties in these numbers are referred to the low^er phase. Here 
a is the negative thermal expansion coefficient. 3 is the temperature gradient in the 
static state, u is the kinematic viscosit>'. Ap (i.e. — p^) is assumed to be positive, 

thereby excluding the Rayleigh Taylor instability, g is the gravitational constant. (Tq 
and Gi are the interfacial tension and its temperature coefficient, respectively. 

The governing equations are nonlinear and admit the conductive (quiescent state as 
a trivial solution. W'e linearize the equations about the trivial base state, eliminating 
all of the dependent variables in favour of the vertical component of velocity and 
temperature. Linearization of a dependent variable A gives 


A = ( 3 . 2 ) 

where A^ is the quiescent state and r is a deviation from this state. A' is further 
decomposed as 

A^ = Ao{z)c'^^P'^3 (3.3) 

This means that a Fourier transform in the ,r-direction and a Laplace transform in 
time has been used. The dependent ^'ariables are the transformed temperature 6 ^q. 
the transformed vertical component of velocity Wq and the transformed interfacial 
deflection 7 / 0 . We get the momentum and energy equations for the upper phase rus 
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n{D^-J)e,r = gO^-^-]^. ( 3 . 5 ) 

Here D reoreseiits the total derivative with respect to 2 . r, s. a. ri and rn are the 
denshy: k.nenrat.c v.scosity, thernral expan^^^^^^^ 

interface are 

u 0^ = gvo = ^'0 ’ 

[DWn\t = 0, 


(3.6) 

(3.7) 


o^)DWo + 2K[s]u;"T>l'l^|^ = {RC 

+ G + u;2)^77o, 


(3.8) 

)2 = iVW(r?o - Qq")- 

(3.9) 

! and heat flux at the interface are 


|0o - /v[l/H^ol- = 

(3.10) 

[K[m]Deo\t = 0. 

(3.11) 

1 J 


K[m] represents an operator that takes the value oi m m me upH- p— ^ 

unity in the lower phase. pi„pnvalne nroblem and has 10 dimensionless 

The above system represents an eigenvalne prooiem 


The sim of the quantity uno{^>)/Vo^ ‘"7" . ^ U;^ve 

tells us whether we h^ive upflow “ Tud equal to 0o/n 

downflow at a crest. ^mf^d tells us that we have a hot spot at the 

is a temperature perturbation temperature perturbation 

crest if it ,s positive Here ( 6 - ft )i ^ the f .„„„'hodes- or scenarios' 


crest if it is positive. Here tc- - — modes’ or scenarios' 

at the interface. It is clear that there are four „ 2 and 

and these words will be used interchangeably. These are aepiciea g 

SSw’l«.S;io”s"rdte el«rnSion: uy ;)) 


these Bow scenarios and the scqnenee oi 'y" ' “Jjyy.,,. eicenfunctions ny ( 
level, the total depth d, and the depth ratio ( (= d*/d h The e^ento^^^^ ^ e^t 

The discussion of the physics assumes that q 0 - 

applying a spectral-r method, it wiui shown that Im(q) - 0 ^^hen He(q) 

cases studied. t pet pH hv recovering the results of 

The reliability of the numerical procedure \\as tested b\ & 
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Figure 2. Four How modes at the fluid -Huid interface. When the interface is flat, there is no 
difference between mode I and IF or mode III and W. 


Table 1. Ph ysical properties of fluids 





air 

above 
water or 
D.C. oil 

water 




L/ow v^oiumg on 

under 

benzene 

(16^0 

under 

air 

(0+°C) 

benzene 

above 

water 

gallium 
under 
D.C. oil 

parameter 

(units) 

under 

air 

above 

gallium 

density 

(g 

0.968 

0.968 

0.0012 

0.998 

0.999 

0.884 

6.09 

negative thermal 
expansion coefficient 

(xio^-* ’r) 

9.6 

9.6 

54.0 

2.06 

-0.68 

14.5 

1.0 

thermal conductivity 
(xlO^ erg enr^ s“^ 

1.55 

1.55 

0.26 

5.97 

5.97 

1.64 

334.0 

thermal diffusivity 
(xlO”'^ cm“ s“^) 

1.1 

l.i 

160 

1.43 

1.43 

1.04 

146.0 

kinematic viscosity 
(stoke) 

1.0 

0.05 

0.152 

0.01 

0.01 

0.0067 

0.003 54 

interfacial tension 
(dyne cm”^ ) 

20.9 





32.8 

74.9 



718.0 

negative inter facial 
tension gradient 
(xlO“^ dyne cm-1 °C) 

5.8 

— 

— 

5.66 

14.0 


38.89 


Zeren Reynolds (1972) and Perm Sz Wollkind (1982). Table 1 contains the fluid 
properties used in our calculations. The results for the silicone oil-air system agree 
with those of Perm & Wollkind (1982) (witliin the round-off error). The comparison 
for the water- l)enzene bilayer shows good agreement with the values of Zeren & 
Reynolds (1972). We note that Renardy Joseph (1985) and \\5ihal &z Bose (1988) 
have independently verified these latter results. 
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4. Discussion of the numerical results 

As we have used a linearized model, the calculations can only give information 
about the flow state at the onset of convection. The discussion is divided into several 
parts, each referring to the physical situation on hand. These are: (a) a liquid-gas 
case, as exemplifled by the silicone oil-air system: (6) a liquid-liquid case, as depicted 
by the water-benzene system: (c) a case of two liquids in the presence of solidification, 
which is depicted by the gallium-silicone oil system with a solidifying gallium phase 
below; and [d) the situation that arises when the temperature coefficient of density 
is positive, as exemplified in the w'ater air system with the temperature range of 
water between 0 and 4°C. As noted by Chandra & Holland (1983), there are some 
commercially important liquid semiconductors, such as mercury cadmium telluride, 
with positive temperature coefficients of density. It is clear from the choice of systems, 
where the upper layer is less dense than the lower layer, that the Rayleigh Taylor 
instability will be excluded from this study. Moreover, as neither phase is in motion 
in the base state, the Kelvin-Helmholtz instability is also excluded. 

For a chosen system, the only control parameters in experiments are the total 
depths, depth ratio and gravity level. The critical temperature difference and critical 
wavelength at the onset are results that come from the linearized stability calculation. 
Before we discuss the results any further, we will clarify the roles of gravity, total 
depth and the parameter The lowering of gravity has the role of increasing the 
relative importance of Marangoni to Rayleigh effects and also reduces the Bond 
number or the effect of gravity waves. In this paper, different gravity levels will be 
chosen and a low^ering of gravity will therefore reduce both gravity waves as well 
as the Rayleigh effect. Unlike Smith (1966), we will not study the case where only 
capillary and gravity waves are considered and where buoyancy is ignored. If the 
total depth is reduced for a given system, we have the effect once again of increasing 
the Marangoni effect relative to the Rayleigh effect in both phases. However, for a 
fixed total depth, decreasing ^ has the role of increasing the Rayleigh effect in the 
lower phase at the expense of the Rayleigh effect in the upper phase. Besides, it also 
has the effect of increasing viscous resistance in the upper phase and this will play a 
role in the flow structure that the fluid bilayer system settles into. In what follows, 
w^e shall refer to figure 2, which shows four possible flow 'modes' or ‘scenarios* at the 
interface at the onset of convection. 

The actual ‘mode' that a s}''stem settles into depends upon the thermophysical 
properties, gravity, total depth and depth ratio, but we will be less concerned about 
the mode that is realized and more concerned with the sequence of transitions from 
one mode to another as control parameters change and. therefore, will make general 
statements regarding this transition sequence for a variety of configurations. These 
statements wall then be 'verified* by appealing to specific calculations on particular 
fluid -flu id systems and do not entail a mathematical proof. 

It is instructive to note the relation betw^een these modes. When the onset state of 
motion of the fluid bilayer goes from mode I to mode II. it simply means that a crest 
must gradually transform into a trough and the flow' directions do not change at any 
particular lateral position along the interface. In going from mode II to mode HI, two 
things happen: first, the crests and the troughs become progressively smaller and, 
second, the interface approaches the base state temperature. This change continues 
until the positions w'hich were formerly crests now' become troughs and vice versa. 
Meanw'hile, the temperature perturbations also reverse in sign (cf. figure 2). Mode IV 
is seen only in one of the cases that w^e discuss, and changes to mode I under some 
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conditions. The transition between I and IV is similar to that between II and III. It 
is clear that the mode transitions must be smooth in the sense that the only way 
in which any two modes can coexist is in the asymptotic case where the interface is 
flat. Otherwise sudden mode transitions would imply the existence of codimension 2 
points. It is easily seen that such codimension 2 points are precluded as the critical 
Marangoni number is simple, unique and obtained as the ratio of two determinants 
resulting from the imposition of the boundary conditions, while the corresponding 
eigenfunction is unique and this issue is made clear by Nadarajah & Narayanan 
(1987). In other words, we cannot ever get two or more coexisting flow modes at the 
critical Marangoni number for the laterally unbounded case. 

There are some specific characteristics of the Marangoni and Rayleigh effects which 
are worth stating. First, if the Marangoni effect is operative alone then hot fluid at 
the interface must move towards the cold spot and the flow mode may be either I 
or II. Now, whether the hot region at the interface is a trough, as in mode I, or a 
crest, as in mode II, depends on the magnitude of the forces and the mechanical 
and viscous resistances in both phases. This is so as the interface at the hot spot 
will then bump towards the region that exerts the greater resistance to flow in that 
region. Second, if the Rayleigh effect is operative alone then matters become a little 
complicated. In the case of a liquid-liquid bilayer, the upper fluid offers resistance 
on account of its viscosity and density and yet it conducts heat. We can expect to 
see any of the modes depending on the forces at play and the magnitude of the 
resistance. As we continue to consider the pure Rayleigh effect, but now restrict the 
study to the 'heated from above' problem, it initially appears that no steady flow 
will occur unless we have the odd case of a positive thermal expansion coefficient. 
However, this premise can be shown to be false. We refer the reader to Gershuni 
& Zhukovitskii (1980). These authors discuss a case where the upper fluid is less 
conductive than the lower fluid and has also a much smaller thermal expansivity. 
In that peculiar case, we can see that a mechanical perturbation to the upper fluid 
sends a hot fluid element towards the interface, where it easily transmits heat to 
the more conductive lower region. This in turn excites buoyancy- driven motion in 
the lower phase and, as a result, momentum from the lower phase is transmitted to 
the viscous upper layer and the process continues. Because Gershuni & Zhukovitskii 
(1980) studied the case of water and mercury, we verified their results as a test of our 
numerical method, but otherwise did not consider this particular system further. In 
the unusual case where the lower fluid has a positive thermal expansion coefficient, 
we obtain mode IV for a liquid-gas system because it is easier to push light cold 
fluid upwards and towards the interface. We will now discuss these problems in the 
following sections and provide numerical evidence for various system calculations. In 
what follows, both buoyancy and interfacial tension gradients, in general, come into 
play unless noted otherwise. 


(a) Liquid-gas system. 

It is clear in the liquid-gas case that only modes I and II are possible candidates, 
as the upper fluid is virtually passiv^e, offers little fluid mechanical resistance and 
hot fluid from below' must flow upw^ards for the reason that this is the situation 
favoured by both buoyancy and interfacial tension gradient forces. Calculations that 
treat the upper gas as passive w^ere compared with those that treat it as active and 
we obtained results that were within 1% of each other insofar as the values of the 
critical temperature difference w^ere concerned. 
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Figure 4. Mode switching for the liquid -liquid (water beiizene) system 'heated from below': 
(a) the effect of gravity level and total depth {(. = 4); (b) the effect of depth ratio (ly). 

in figure 2, and in the liquid-gas system we can only obtain the first two modes. 
Ferm & Wollkind (1982) indicated that the mechanism changed from Marangoni- 
to Rayleigh-dominant when there was a change in the slope of the graph between 
critical temperature gradient and lower phase depth. We feel that our classification is 
more definitive as it identifies various scenarios at the interface between both fluids. 

(b) Liqmd-liquid systems 

Liquid -liquid systems are interesting as they may be examined with two different 
heating directions and we choose the water- benzene system to explain the physics 
because it provides a good test case to verify the detailed results of Zeren & Reynolds 
(1982). 

The modal sequence as we increase gravity, and also ‘heat from below', goes from 
mode I to mode II and then to mode III. The reason why we can expect to see 
mode III as gravity is increased from mode II is because the Marangoni effect becomes 
less important and less work is required to push up the short cold columns seen in 
III as opposed to the tall cold columns seen in IV. In fact, mode IV occurs for the 
peculiar situation where the thermal expansion coefficient is positive, and this we will 
discuss later. The particular system that we have investigated has a large value of a, 
the thermal expansion coefficient ratiof, and therefore it is generally biased towards 
greater buoyancy in the upper phase in comparison to the lower phase, unless of 
course i becomes so small so as to discourage buoyancy in the upper phase in favour 
of the lower phase. It is for this reason that, as gravity is increased from mode II, 
the fluid goes into mode III for all of the i values that we have used. Further, it is 
at the point of mode switching that the Marangoni effect ceases to be of importance 
and this is seen by the fact that the fluid at the interface does not move from hot 
regions to cold regions. However, once the fluid is in mode III. the Marangoni effect 
delays the instability by raising the critical temperature difference. When the mode 
changes from mode II, upon decreasing gravity, then mode I. which is favoured by 
the Marangoni effect, is realized. It should be noted that the mode can remain in 
mode II and never change into mode I on the reduction of the gravity level. This 
peculiarity once again has to do with the value of i and upon the viscosity ratio 


t All of our calculations and statements in the Ccxse of liquid -liquid systems are restricted to the 
situation where the ratio of the upper phase thermal expansion coefficient to the lower one is much 
greater then unity. The reverse case is not discussed in the cause of brevity but the physically based 
arguments follow in a similar manner. 
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s. If the upper fluid is very viscous and / is small then more resistance to flow is 
exerted by the upper phase, preventing the interfacual tension from flattening out a 
crest and thereby disallowing a mode II structure to de\^lop into a mode I type. It 
turns out that this observation was mathematically proven by Smith (1966) for the 
problem that considered only gravity and capillary waves without buoyancy effects. 
We have also checked and verified Smith s assertion by considering a water benzene 
bilayer with a value of f — 0.5 for very low gravity levels and obtained a mode II 
flow structure. 

The roles of total depth and f are somewhat complicated for the liquid-liquid 
problem that is ‘heated from below'. As tlie effect of reducing the total depth is 
to favour the Marangoni over the Rayleigh effect, this simply means that the mode 
switching for this heating direction goes from mode III to mode II and then to mode I 
as total depth is reduced. Figure 4n bears this out. Of course, for very low T mode II 
is obtained on account of the viscous-mechanical resistance to flow in the upper phase 
and figure 4b bears this out. 

Now. if total depth and gravity are kept constant and ( alone is increased, physical 
reasoning demands that the buoyancy effect in the upper layer increases relative to 
the lower layer. Let us suppose that the fluid properties and conditions are such that 
the fluid settles into mode II. If the total depth is smallf and we continue to increase 
L we expect the flow to switch from mode II to mode I, because the interface gets 
closer to the hot lower surface in what is already a thin layer and the Marangoni 
effect plays a dominant role when the interface gets close to the lower hot surface, 
thereby giving rise to a hot trough. Meanwhile, in this thin Iciyer. the large value of f 
causes the Rayleigh effect in the upper layer to become more significant than in the 
lower layer. Buoyancy in the upper layer tends to cause hot plumes to rise in that 
layer, in opposition to the Marangoni-iiifluenced flow at the interface. As a result, we 
see a vertical stacking of flow cells in tlie upper phase as f is increased. The upper 
cell in the upper phase is a result of the buoyancy in tliat j)hase and the lower cell 
in the upper j)hase is a result of the Marangoni motion tliat results on account of 
the proximity of the interface to the lower heated surface. The cell stacking persists 
even as the mode switches from II to I. The effect of / for a small total depth is seen 
in the left region of figure 45. while figure 5 depicts the velocity and temperature 
eigenfunctions when vertical stacking takes place. 

Contrast the situation discussed above with the case where the total depth is large: 
we then expect the fluid to switch from mode II to mode III as / is increased and 
the explanation for this is as follows. The large total depth increases buoyancy in 
both phases; the upper phase being even more buoyant than the lower on account 
of two reasons. First, the increasing value of f continues to enhance the buoyancy 
in the upper phase in relation to the lower and second, but le.ss important, the 
value of a is much greater than unity (in the systems that we chose to compute). 
As the upper layer is now very buoyant, it would require less work for the fluid 
to go into mode III, where a shorter column of fluid is pushed up against gravity 
in both phases, in contrast to the situation of mode II. Remember that a shorter 
column of liquid is puslu'd up against gravity in mode I as well, but the system 
would prefer mode III as the upper phase portrays more buoyancy in this mode as 


t The notion of large total depth or small total depth is entirely system dependent and ordy calcula- 
tions or experiments can determine how large the total depth should be for us to s<»e what we predict. 
Eiowever. we can still make statements of a qualitative manner and verify them by calculations. 
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Figure 6. Mode switching for the liquid-liquid (water-benzene) system ‘heated from above’: 
(a) the effect of gravity level and total depth (dt = 2 mm): (6) the effect of depth ratio (1^). 


compared to mode L The numerical calculations confirm the physical arguments, as 
seen in the right-hand region of figure 4b. In these cases, where a large total depth is 
considered, an increasing value of £ does cause the interface to get closer to the hot 
lower rigid surface, however it is not close enough to encourage an overriding effect of 
Mar angoni- influenced motion, and that is why the system settles into mode III with 
cold troughs. At the intermediate values of total depth, mode II remains as expected 
for all values of 

When the water-benzene bilayer is Tieated from above*, we predict a sequence 
from mode I to mode II as we increase gravity; the explanation for this mode switch 
is as follows. ‘Heating from above' causes convection that is started by a Marangoni 
influence and therefore the flow must necessarily be in either mode I or mode II. 
Meanwhile, gravity serves the purpose of stabilization and therefore delays the in- 
stability. Now the value of a, the thermal expansion coefficient ratio, is much greater 
than unity in the system studied and an increase in gravity causes an increase in the 
resistance to flow in the upper phase because of the stabilization effect. This in turn 
causes the switch into mode II, creating hot crests at the interface. 

Likewise, for the ‘heated from above* case, the mode switching goes from mode II 
to mode I as we decrease total depth for the reason that a decrease in total depth 
is tantamount to a decrease in the stabilizing Rayleigh effect in each phase, thereby 
enhancing the Marangoni effect. If we again consider gravity and total depth to 
remain constant but vary £ for the ‘heated from above* case, we will go from mode II 
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to mode I as ^ is increased. This mode switch will take place when the total depth is 
small, because an increase in ( has the effect of increasing the viscous resistance in the 
lower phase, even though it makes the upper phase more stabilizing than the bottom 
phase. When the total depth is large, the effect is mainly to make the upper phase 
stabilizing, thereby increasing the resistance there and so mode II remains intact for 
all £. Note, as a result of our reasoning, the ‘heated from above' configuration cannot 
cause vertical stacking. Vertical stacking can occur only when both Marangoni and 
Rayleigh effects destabilize in a particular phase and yet when both act in opposition 
to each other in a 'flow direction’ sense. This does not occur in the 'heated from 
above’ configuration with a negative thermal expansion coefficient. Our thinking 
in the above is underscored by the numerical calculations that are presented in 
figures 6a, 6. 

A comment on Smith’s (1966) paper is in order. His study was concerned with 
Marangoni convection in the presence of gravity waves but without a buoyancy mech- 
anism. No mention was made of hot or cold spots and only modes where fluid from 
the lower phase moved upwards or downwards from a crest were considered. Yet a 
pair of sufficient conditions were obtained to indicate whether flow moved upwards 
or downwards into a crest. This was obtained by use of a flow indicator as in this 
study. All of our observations validated Smith’s (1966) derived result. 

(c) Solidifying phase below a liquid-liquid bilayer 

The third case is a modification of the 'liquid-liquid’ problem that is 'heated from 
above’. It involves a solidifying phase below the lower layer of liquid. The condition 
at the boundary of the lower liquid and the adjacent solidifying phase is replaced by 

DOq = coth(a;yl). (4.1) 

It is implicitly assumed that the rate of solidification is slow enough to be negligible 
and so that there is no net flow in the base state. The perturbed deflection of the 
solid-liquid interface is given by 

Co = Oo (4-2) 

In the above, A is the dimensionless thickness of the lower solidifying phase, scaled 
with respect to the lower liquid thickness. We note that our calculations exclude the 
important situation where constitutional supercooling, as considered by Mullins & 
Sekerka (1964), is involved. The main result from our calculation is that the solid 
thickness of the lower phase does not affect the flow structure at Ig but does affect 
the flow at low gravity levels. In other words, the coupling is onh^ one way at high 
gravity. We surmise that even though the solid thickness destabilizes the flow, it is 
insignificant compared to the overwhelming stabilization of gravity. This result is 
depicted in figure 7. It is apparent from this figure that solidification at lower levels 
of gravity may not prove beneficial at all under this heating arrangement as the 
critical temperature gradient is lowered in the low gravity state and convection is 
thereby enhanced. As a result, it appears that liquid-encapsulated crystal growth 
is better conducted under Earth’s gravity conditions as the configuration would be 
more stable. 

We can also determine whether the solid-liquid interface deflection is in phase or 
out of phase with the upper interface by examining the ratio C( 0 )/^( 0 )- What is even 
more interesting is that when mode switching takes place as we increase gravity, we 
find that the deflection at the solid-liquid interface lower phase does not change sign. 
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Figure 8. Mode switching for the liquid -gas (water-air) system 'heated from above’. The liquid 
in this system has a positive thermal expansion coefficient (water is between 0 and 4°C): (a) 
the effect of gravity level and total depth {( = 1); (b) the effect of depth ratio (Ig). 
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{d) The case of a system with a positive thermal expansion coefficient 

Finally, we consider the case of a system where the thermal expansion coefficient 
of the lower phase is positive. This is exemplified by the results using the water air 
system, as shown in figures 8a, b. Here the heating can take place from above or 
below. If the heating is from above, the origin of the convection can only be due 


- I I i I 

I I I 1 


Proc. R. Soc. Land. A (1995) 



Bilayer Rayleigh- Marangoni convection 
Table 2. Summary of mode switching 


501 


liquid-gas system 
heated from below' 


negative thermal expansion 
coefficient in liquid 

liquid-liquid system liquid-liquid system 

‘heated from below' ‘heated from above’ 


increasing 

gravity level I — *■ II 

increasing 

total depth I — ^ II 

increasing 

depth ratio £ II — I 


I ^ II III 
I II ^ III 

small total depth large total depth 
II — I II ^ III 


I II 

I ^ II 

II I 


posit iv^e thermal expansion 
coefficient in liquid 


liquid-gas system liquid-gas system 
‘heated from above' ‘heated from below’ 


increasing 




gravity level 

I- 

^ IV 

I 

increasing 
total depth 

I - 

- IV 

I 

increasing 




depth ratio i 

IV 

— I 

I 


to buoyancy as the upper layer is virtually passive. The flow, as expected, stays in 
mode IV for large gravity lev^els. While the convection is necessarily of buoyancy 
in origin, it is also true that the Marangoni effect plays a part in countering the 
flow, thereby delaying the instability. Since the strength of this interfacial mode of 
convection is dependent on the magnitude of the interfacial tension gradient, we 
can observe that at low gravity the Marangoni effect is dominant and, as gravity 
or total depths are reduced, the mode switches to I. Here Marangoni convection 
causes fluid to move from hot to cold regions at the interface and a shorter column 
of hot heavier water is balanced by a taller column of cold lighter water. As gravity 
decreases, the system becomes more stable (despite the change in the flow mode 
from IV to I) as evidenced by an increase in the critical temperature gradient, just 
as we would expect. If the water-air system were ‘heated from below’ so as to allow 
the temperature range across the water to be between 4 and 0°C, then it is likewise 
argued that we would get mode I at all gravity levels, as shown in figure 9. Here the 
buoyancy stabilizes and the onset flow comes from the interfacial tension gradient. 
Much like the ‘heated from below* configuration in the liquid-gas case for negative 
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thermal expansion coefficient and pure Marangoni flow, mode I will be the structure 
choscai by the system. 

Tai)ie 2 gives a summary of the mode switching explained in the study. Once again, 
our results depend ui)on the thermophysical properties of the individual systems 
chosen in this study. However, the physically based arguments can be similarly made 
for otlu^r systems. We believe that it is valuable to make such arguments in order to 
identify the mechanisms that dominate the flow structure at an interface. 
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Marangoni convection in multiple bounded 
fluid layers and its application 
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A brief review of multilayer convective phenomena that is associated with materials 
processing is presented. Several instability phenomena that can occur in a bilayer of 
two fluids heated from either above or below and the effect of laterally and vertically 
confined geometries are explained. In particular it is shown that such confinement 
can lead to the occurrence of codimension-two points and pure thermal coupling that 
is initiated by convection in an upper gas phase during liquid-gas bilayer convection. 
Experimental evidence that shows the effect of geometrical restrictions is given. 

Keywords: Rayleigh; Marangoni; interfacial tension driven convection; 
buoyancy driven convection; multiple fluid level convection 


1. Introduction and physics 

Much of the work reported in this paper has been motivated by the need to under- 
stand a technique for growing certain crystalline materials, known as the liquid- 
encapsulated vertical Bridgman (LEVB) crystal growth method. Liquid-encapsulated 
crystal growth is a process for producing III-V semiconductor crystals from bulk liq- 
uid melts. The demand for crystals of increasingly higher purity and lower defects 
requires us to understand this process in much greater detail. Some examples of crys- 
tals grown using this technique are gallium arsenide (GaAs) and indium phosphide 
(InP). Taking GaAs as an example, when GaAs is melted, it has a tendency to decom- 
pose, releasing arsenic gas and destroying the desired stoichiometric ratio. To prevent 
this decomposition, a liquid encapsulant of boric oxide (B2O3) is placed on top of 
the gallium arsenide. In addition, an inert gas may be placed on top of the B2O3. 
These three layers are placed in a crucible, w^hich is low^ered through a temperature 
gradient created by a furnace. The lower end of the crucible is cooled, thereby solidi- 
fying the gallium arsenide. This configuration is shown schematically in figure 1. The 
heating configuration generates vertical as well as radial temperature gradients and, 
consequently, interfacial-tension-gradient-driven convection, also known as Marango- 
ni convection, and buoyancy-driven convection, also called Rayleigh convection, can 
occur at the liquid-gas and liquid-liquid interfaces as well as in the bulk fluid regions. 

While the LEVB technique is the motivation for this study, only by considering 
simple systems can we have a clearer understanding of the physics of the convective 
process. Radial gradients of temperature, creeping of the encapsulant along the ver- 
tical sidewalls and solutal gradients all have a complicated effect on the convection. 
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Figure 1 . Schematic of a liquid encapsulated crystal grower: a system of three convecting fluid 
layers. Convection in the GaAs liquid influences the quality of the GaAs solid. 

Indeed tlie onset of convection in an actual LEVB system occurs simultaneously 
with the application of any temperature gradient. However, a clear understanding 
of comection in LEVB and many other materials processing methods requires us 
to consider problems where classical fluid mechanical procedures may be employed, 
thereby- simplifying the mathematics while simultaneously revealing the essential 
physical h'atures. 

One such problem is the Rayleigh-Marangoni problem. Here a gas or another liquid 
superimposes a liquid layer and a vertical temperature gradient is applied. Suppose 
that the density and interfacial tension of the liquid decreases with increasing tern** 
perature. As the liquid is heated from below, it is top heavy. A small disturbance 
can upset this arrangement if the overall temperature difference is large enough and 
flow can ensue in the form of buoyancy or Rayleigh convection. However, flow can 
occur even in the absence of gravity. For example, in the quiescent state the liquid- 
gas interface is flat and a small disturbance to it causes a transverse temperature 
gradient at the interface causing fluid to flow from warm regions of low interfacial 
tension to cold regions of high interfacial tension. Hot fluid from below rises to the 
interface and cold fluid from the interface moves down to maintain continuity of fluid 
flow and the convection continues as Marangoni convection. For small values of the 
vertical temperature gradient, the fluids remain quiescent and transport heat by pure 
conduction. However, when the temperature gradient reaches a critical value, even 
the smallest disturbances imposed on the system amplify with time and the system 
reaches a steady or time-periodic steady state. In other words a critical temperature 
gradient is required for convection to occur. More details on the nature of this type 
of convection are available in the reviews of Koschmieder (1993) and Davis (1987). 
We will explain the physics of single and multilayer convection in laterally bounded 
geometries where the layers are heated from above making them gravitationally sta- 
ble and where the layers are heated from below making them gravitationally unstable. 
The explanation of physics in multilayers will be followed by a report on two sets of 
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Aspect Ratio 

Figure 2. Plot of the critical Marangoni number versus the aspect ratio of a cylinder. The mode, 
m, with the smallest Marangoni number at a given aspect ratio is the mode or flow pattern at 
the onset of convection. 

experiments. Noticeably absent from this paper will be the effects of solutal convec- 
tion, some aspects of which have been covered by several other authors (McFadden 
et al. 1984; Turner 1985). It should be noted that this paper is a brief report of the 
work done by us and a few other researchers. Greater detail is available in the thesis 
by Johnson (1997) as well as papers by Johnson & Narayanan (1996, 1997, 1998). 
Fuller explanations of the effects of convection on crystal growth are given by Hurle 
(1994), Muller (1988) and Schwabe (1981). 

The extent of convection is often characterized by a dimensionless temperature 
difference represented by the Marangoni or Rayleigh numbers. The Marangoni num- 
ber is proportional to the depth of the liquid, the temperature difference and the 
variation of the surface tension with respect to the temperature and inversely pro- 
portional to the dynamic viscosity and thermal diffusivity. The Rayleigh number 
is proportional to the cube of the liquid depth, gravity, the temperature difference 
and the thermal expansion coefficient and inversely proportional to the kinematic 
viscosity and the thermal diffusivity. In a physical system, fixing the temperature 
difference necessarily fixes both the Rayleigh and Marangoni numbers. 

We begin by confining our discussion to a single layer of fluid, heated from below, 
with a free surface. In this configuration both buoyancy and interfacial-tension forces 
become important. For larger depths, buoyancy is more important than interfacial- 
tension effects, and when the fluid depth is small, interfacial-tension forces dominate 
convection. 


(a) Physical effects of a bounded geometry 
Consider a single layer of fluid bounded below by a rigid conducting plate and 
whose upper surface is bounded by a passive gas. By a passive gas we mean a. gas 
which has no viscosity and only conducts heat away. The lower plate here is at 
a higher temperature than the passive gas. In a fluid of infinite horizontal extent. 
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Figure 3. Schematic of four different flow patterns: o. fluid flowing up; x, fluid flowing down. 


there is no limit on the number of convection cells. However, in a bounded finite-sized 
container only a finite number of convection cells may exist. Physically this means 
that at the onset of convection in a bounded cylinder, only one flow pattern will 
usually exist. As the aspect ratio (radius/height) of the container increases, more 
convection cells will appear. Figure 2 is a representative calculation of the critical 
Marangoni number for various aspect ratios and for different azimuthal modes m. 
The Biot number, which is a dimensionless surface heat-transfer coefficient is equal 
to 0.3. 

In a bounded cylinder, each flow pattern is associated with an azimuthal mode, 
m, and radial mode, n. For example, at an aspect ratio of 1.0 in figure 2, there is 
an m — 1. u — 1 flow pattern (see figure 3). For an aspect ratio of 1.5, there exists 
an m = 0. n = 1 flow pattern. The azimuthal mode is the number of times the 
azimuthal component of velocity goes to zero, and the radial mode is the number 
of times the radial component of velocity goes to zero starting from the centre for a 
given vertical cross-section. 

At particular aspect ratios, where the fluid switches from one flow pattern to the 
next, there coexist two different flow patterns. These aspect ratios are known as 
codimension-two points. For certain codimension-two points, the flow patterns will 
interact nonlinearly to yield oscillatory behaviour (Rosenblat et al 1982; Johnson & 
Narayanan 1996). This phenomenon will be shown later in §2. 


(b) Physical effects of multiple fluid layers 

Imagine a less dense immiscible layer of fluid above the lower layer of fluid. Here the 
lower layer is bounded below by a rigid conducting plate and another rigid conducting 
plate bounds the upper layer. Once again let the temperature of the lower plate be 
greater than the upper plate. The interface between the two fluids may deform and 
is capable of transporting heat and momentum from one layer to the other. We will 
now consider the various types of convection that can occur in a bilayer of two fluids. 

In order to distinguish the various convection mechanisms, w^e introduce phrases 
such as ‘convection initiating in one layer or another'. Strictly speaking, convection 
occurs in both fluids simultaneously, although one layer may be more unstable than 
the other, driving flow in the other layer. 

Turning now to various convective mechanisms, consider figure 4. Suppose that 
convection initiates in the lower layer. The upper layer responds by being dragged. 


Phil. Trans. R. Soc. bond. A (1998) 



889 


Marangoni convection in multiple bounded fluid layers 


{a) ih) 


GO GO 

0000 

OOCIO 

0000 


ic) id) ie) 


0000 

00 

[JO 

Ot 

300 

0000 

oooo 

oooo 


Figure 4. Schematic of the different types of convection-coupling: (a) lower dragging mode; (b) 
viscous coupling; (c) thermal coupling: (d) upper dragging mode; (e) pure thermal coupling. 
Moving from (a) to (e), the buoyancy force in the upper layer increases. Gas-liquid thermal 
coupling, with surface-driven flow, is caused by the upper fluid buoyantly convecting and simul- 
taneously inducing interfacial-tension- or buoyancy-driven convection in the lower layer near the 
interface. 


generating counter rolls at the interface. Hot fluid flows up in the lower layer and 
down in the upper layer. The upper layer is not buoyant enough and moves by a 
combination of viscous drag and the Marangoni effect. This is seen in figure 4a. 
The sign of the velocity switches and the maximum absolute value of the lower-layer 
velocity is much greater than the maximum absolute value of the velocity of the 
upper layer. 

When the buoyancy in the upper layer increases and the upper layer begins to 
convect, one of two things can happen. The first possibility is that the two fluids are 
inscously coupled. Physically this can be shown in figure 46 as counter-rotating rolls 
in the two fluids. This can also be denoted by the vertical component of velocity 
switching sign at or near the interface, while the temperature perturbations indeed 
switch sign at the interface itself. If the temperature perturbation switches sign 
near the interface in either layer near the interface we would say that the bilayer is 
nearly viscously coupled. In particular if the switch takes place in the upper fluid 
near the interface, then the lower layer is slightly more buoyant. If the temperature 
perturbation switches sign in the lower layer, then the upper layer is more buoyant. 
The Marangoni phenomenon, for fluids, whose inter facial tension decreases with an 
increase in temperature, plays an ambiguous role here. The hot fluid flowing up in the 
lower layer causes the fluid at the interface to move in the same direction. However, 
the colder fluid moving down in the upper layer contradicts this. The exact effect 
the Marangoni phenomenon has on the two fluids depends on where the thermal 
perturbations change sign. For the situation where the thermal perturbations switch 
sign at the interface there is no Marangoni effect. 

The second possibility is thermal coupling where the rolls are corotating (figure 4c). 
Here hot rising fluid from the lower layer causes hot fluid in the upper layer to 
flow up. The maximums of the vertical component of velocity and the temperature 
perturbations have the same sign in each fluid layer. Strictly speaking, the transverse 
components of velocity should be zero at the interface. However, thermal coupling is 
sometimes referred to the case when a small roll develops in one of the layers so as 
to satisfy the no-slip condition at the interface. In this situation, when the interfacial 
tension decreases with an increase in temperature, the Marangoni effect encourages 
flow in the lower fluid layer, and discourages the flow in the upper fluid. 

Another interesting phenomenon is present at certain fluid depths where both ther- 
mal and viscous coupling can occur. At these depths, a competition arises between 
the two types of convection. As both convection configurations cannot occur simulta- 
neously, the fluids begin to oscillate between these two states. This phenomenon was 
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Figure The four possible iiiterfacial structures at a fluid fluid interface. Each structure can 
give information about the driving force of the convection. 


first reported by Gershuni & Zhukhovitskii (1982). and has recently been confirmed 

bv Anclereck et ai (1996). , • -j.- 4 . 

As tlu' buovanev continues to increase in the upper layer, convection initiates 

in onlv the upper layer and the lower layer is viscously dragged (figure 4d). This 
situation onlv occurs when the upper huid is a liquid, as gases are very tenuous and 
will not ('xert much shear. The vertical component of velocity m this case switches 
sign and the magnitude of convection in the upper fluid is much greater than the 

magnitude of convection in the lower fluid. 

The last figure (figure 4e) is an example of what may be called pure thermal 
coupling. This typically occurs in a liquid-gas system where iiuoyancy convection 
is predominant in the gas layer. The convecting gas then simultaneously creates a 
non-uniform temperature profile across the liquid-gas interface and generates either 
Maraugoiii or buoyancy-driven convection in the lower layer (. Johnson et af. 1998) 
Notice that the convection in the lower layer is now generated purely by horizontal 
temperature gradients at the interface and not by viscous dragging. To maintain the 
no-slip (‘oudition at the interface a small counter-roll may develop m the gas-phase. 
This roll is not shown in figure 4e. 


(c) Physics of interfacial structures 

Another indicator of what is occurring in bilayer convection can be inferred from 
the fluid-fluid interface instead of the bulk convection. In a paper by Zhao et al. 
(1995) four different interfacial structures were identified for any given convecting 
bilaver with a deflecting interface. Each of these structures depends upon whether 
fluid was flowing into or awav from the trough or the crest, and whether the fluid 
was hotter or cooler at the trough or the crest of the interface. Hot fluid flowing into 
a trough defines the first interfacial structure. The second mterfacial structure has 
hot fluid flowing into a crest. The third structure has hot fluid flowing aw^ froni 
a crest and the fourth structure has hot fluid flowing away from a trough. Each of 

these four scenarios is given in figure 5. . i i* x* r 

One of the important factors to consider in interfacial structures is the direction ol 

the flow along the interface. As interfacial tension is usually inversely proportional to 
temperature, at cooler regions of the interface, the interfacial tension will be higher 
and will pull on the interface. Where the interface is hotter, the mterfacial tension 
will be lower causing the fluid to move away from warmer regions. Another important 
factor is the direction of the flow into or away from a crest or a trough. One re^on 
the interface deflects is due to bulk convection, caused by buoyancy effects, pushing 
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against the interface. Consider two fluids wliose dynamic viscosities are equal. If 
buoyancy-driven convection is occurring mostly in the lower layer, then the fluid will 
flow up from the lower layer into a crest. If the fluid flows down from the top layer 
into a trough, then one would argue that buoyancy-driven convection occurs mostly 

in the upper fluid. . 

In each of the four cases, the interfacial structure can be used to indicate the 
driving force of the convection. In the first interfacial structure, the dominating 
driving force is interfacial-tension-gradient-driven convection. This is seen as the cold 
fluid, with the higher interfacial tension pulling the fluid up into the crest. The first 
inter’facial structure can also occur by buoyancy-driven convection in the upper layer, 
when the density of the upper layer increases with an increase in temperature. In 
the second interfacial structure, buoyancy drives convection in the lower phase. The 
hot rising fluid pushes the interface upwards. As the fluid moves along the interface, 
it cools and eventually sinks back down. The third interfacial structure is dominated 
by buoyancy-driven convection in the upper phase or by interfacial-tension-driven 
convection where the interfacial tension increases with respect to temperature. The 
fourth interfacial structure only occurs when the lower fluid has a positive thermal- 
expansion coefficient. In other words, the density increases with an increase in the 
temperature, causing the cooler lower fluid to flow up into a crest. 

Knowledge of interfacial structures will be beneficial in the understanding of cer- 
tain materials processing problems such as drying of films, coatings and deposition. 

( d ) Physics of heating from above 

In the previous subsections we talked about some of the phenomena that occur in 
single- and multiple-fluid layers heated from below. However, in an attempt to avoid 
convection in crystal growth, the crucible is often cooled from below and heated from 
above. This heating configuration changes the physics of the problem, which is the 

topic of this subsection. , , ■ 

When a layer of fluid is being heated from above, it creates a stable density strat- 
ification. Therefore not only does the buoyancy force not cause convection, it acts 
to inhibit other instabilities. Marangoni convection, though, may still occur in fluids 
being heated from above. 

First we will consider a single layer of fluid superposed by a passive gas. it the 
upper gas is truly passive, then pure Marangoni convection will not occur. For exam- 
ple, suppose some random perturbation causes some part of the surface to become 
warmer than the rest of the surface. The interfacial tension will decrease in this 
region and the tension will pull fluid away from this hot spot. By continuity, fluid 
lying below the hot spot will rise up to replace the displaced fluid. As the lower 
fluid is cooler than the surface this region now cools off and the interfacial tension 
increases, thereby restabilizing the region. However, in real systems, the upper fl^uid 
is never truly passive. Given the same scenario, fluid movement along the interface 
will also drag warmer fluid from above. This warmer upper fluid will further increase 
the temperature in this region, and, depending upon the ratio of thermal-physical 
properties, reinforces the instability. 

By this argument, it appears that an active upper fluid is necessary to have 
Marangoni convection when the system is being heated from above. However, this 
is not the case if the buoyancy effects are included. In Rednikov et al. (1998), it 
was demonstrated, theoretically, that oscillatory onset of convection may occur or 
a single layer of fluid with a purely passive upper gas. The explanation is as follows. 
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If a small volume of fluid is displaced within the bulk of the fluid, the density strat- 
ification acts as a restoring force, causing a dampened oscillation within the fluid. 
These are often referred to as internal waves. The Marangoni force acts similarly, as 
discussed above, also giving dampened oscillations. Apparently when these two forces 
act together they can overshoot one another leading to sustained oscillatory convec- 
tion. Indeed, as was demonstrated in their paper, this only occurs in certain fluids, at 
certain depths, where the buoyancy and interfacial- tension forces are approximately 
equal. 

A completely diflferent type of instability is also possible in two layers of fluids 
being heated from above. Gershuni Zhukhovitskii (1981) first demonstrated this 
phenomenon by analysing two immiscible fluids where the interface between the fluids 
was assumed flat and the Marangoni phenomenon was neglected. They found the 
onset of steady convection when the thermal conductivity and thermal expansivity 
of the lower fluid was much greater than that of the upper fluid. 

The mechanism of this instability is as follows. Suppose an element of fluid in 
the upper layer, near the interface, is displaced towards the lower layer. Because 
the thermal expansion of the upper fluid is so small, this element remains in a 
relatively neutrally buoyant state. Also, as the thermal conductivity of the upper fluid 
is small, it cools very slowly. When the two fluids are heated from above, the displaced 
fluid will be warmer than its surroundings. This element of fluid then heats part of 
the lower fluid near the interface. The lower fluid, with its relatively large thermal 
conductivity and thermal expansivity, quickly heats up and expands horizontally. 
This expansion then causes convection in the lower fluid layer and propagates by 
viscously coupling with the upper fluid layer. The Marangoni phenomenon, if it were 
considered, would act to enhance this instability. 

Another case of interest is convection induced by the Rayleigh -Taylor instability. 
This phenomenon can occur in two immiscible fluid layers being heated either from 
above or below, when the densities of the two fluids are approximately the same and 
the thermal expansivity of the lower fluid is much greater than the thermal expansiv- 
ity of the upper fluid. Upon heating, the density of the lower fluid will decrease and 
become less than the density of the upper fluid. Consequently, the heavier upper fluid 
will begin to sink causing large deformations in the liquid -liquid interface, generat- 
ing the Rayleigh-Taylor instability (Chandrasekhar 1961). This problem has been 
investigated extensively in Renardy i.: Renardy (1985) and Renardy (1996), but is 
of application to materials processing only if the densities of both layers are similar. 

2. Some experimental observations 

The experiments were used to investigate both the oscillatory behaviour near 
codimension-two points and the pure thermal coupling of air with silicone oil (see 
figure 4e), Details on the experimental procedure are available in Johnson (1997) 
and Johnson & Narayanan (1996. 1998). 

(a) Experimental apparatus and procedure 

The experiments consisted of two compartments: one for the lower fluid and one 
for the air. Lucite inserts were used to give a variety of different fluid depths and 
aspect ratios. A copper plate was placed below the liquid insert and the air insert 
was bounded above by a high-thermal-conductivity infrared transparent zinc selenide 
(ZnSe) window. Heating of the copper plate was done by an enclosed stirred water 
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bath that was in turn heated by a hot plate. Xhe top of the ZnSe window was kept 
at a constant temperature by accurately controlling the temperature of the overlying 
air. The overall temperature control w^as kept within ±0.05 °C. 

The flow patterns that developed at the silicone oil^air interface layer were visu- 
alized with an infrared camera. Although other flow visualization techniques could 
have been used, such as shadowgraphy or particle tracing, the IR, camera was chosen 
to prove the viability of its use with opaque materials, such as gallium and gallium 
arsenide. The IR imaging technique is also useful in observing weak thermocapillary 
flow near the surface, whereas shadowgraphy requires some strength in the domain 
flow. 

To guarantee that the flow pattern seen was indeed the flow pattern at the onset 
of convection, the temperature diflFerence applied across the bilayer system was care- 
fully increased. At first a temperature difference was applied that was less than the 
critical temperature difference necessary for the onset of convection. This, and all 
temperature differences, were kept constant for several characteristic time constants; 
ca. 3-4 h. If no flow pattern was seen, the temperature difference was then increased 
by as little as 0.05 °C. This was repeated until the temperature profile at the inter- 
face changed to some distinct pattern, indicating that the fluid had begun to flow. 
At this point, the flow pattern was recorded and the temperature difference noted. 

(6) Experimental observation of codimension-two points 

As was noted in § 1, there exist certain liquid aspect ratios where two different 
flow patterns coexist. For example, in figure 2 at an aspect ratio of 2.3, there exists 
a codimension-two point between the azimuthal modes m — 0 and m = 1. The 
questions we want to answer are: What happens at these aspect ratios? Does one 
flow dominate over the other? Do the different flow patterns coexist as a superposition 
of both states, or do they oscillate and interact between these two states? 

Rosenblat et al. (1982) have performed a weakly nonlinear analysis to investigate 
these questions. They found that all three of these possibilities may occur, depending 
on the Prandtl number, the particular codimension-two point being investigated, and 
on which side of the codimension-two point the aspect ratio lies. To simplify their 
calculations, a vertical and tangential vorticity-free side-wall was assumed. Later a 
more realistic no-slip condition was applied (Zaman & Narayanan 1996; Dauby et al. 
1997), where it was noted that the order of azimuthal modes, as the aspect ratio was 
increased, was different than the vorticity-free condition. These latter calculations 
were done assuming a linearized instability analysis. Therefore, a direct comparison 
of the nonlinear analysis with the experiment is not currently possible. Nonetheless, 
some of the qualitative features should still hold true. 

A series of experiments were performed to first find the codimension-two points 
and then determine the flow patterns at or near the codimension-two point (Johnson 
& Narayanan 1996). A 5.0 mm-high 2.5 aspect-ratio liquid insert was used in con- 
junction with a 11.2 mm air height. Table 1 shows the calculated critical Marangoni 
numbers for four different azimuthal modes for a 2.5 aspect ratio. The table predicts 
that an m = 0 flow pattern should be seen at the onset of convection. However, 
the critical Marangoni numbers for m = 1 and m — 2 are also very close to the 
onset point. Physically, this means that for temperature differences slightly above 
the critical temperature difference, these modes may affect the flow pattern. 

At the onset of convection, a very faint m = 0 double-toroidal-flow pattern could 
be seen. When the temperature difference was increased by just 0.05 °C, the flow 
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Table 1. Critical Marangom number for the first four azimuthal modes for a 2.5 aspect ratio. 


modes 

Marangoni 

number 


0 

69.37 


1 

70.84 


2 

70.41 


3 

72.98 




Figure 6. Infrared images showing the mode-switching behaviour in the paper by Johnson N 
Naravaiian (1996). The experiment used 91 cS silicone oil and a 5.0 mm, 2.5 aspect-ratio insert. 


pattern changed from the double toroid to a dynamic switching between two and 
one flow cells (see figure 6). 

At first, two symmetric cells appeared (figure 6a). Then, one of the cells would 
grow and push the other cell out of the picture, forming a superposition of the 
m = 0 and m = 1 flow pattern (figure 6b), Next one cell would grow (figure 6c), then 
split into two cells, rotated by 90° (figure 6d). This process would then repeat itself 
(figure 6e) arriving back to the original rn = 2 flow pattern (figure 6/). As long as 
the temperature difference was held constant- this dynamic process would continue 
repeating itself approximately every 20 min. 

It is noteworthy that oscillatory convection is of particular importance in crystal 
growth. It has been shown (Hurle 1994) that fluctuating temperatures in the liquid 
melt have a deleterious effect on the crystal quality, leading to a higher dislocation 
density. 

(c) Experimental observations of thermal coupling 

The thermal coupling of air with the lower fluid was originally discovered by a 
series of experiments using the same experimental apparatus (Johnson &: Narayanan 
1998). As was explained in §1. air can thermally couple with the lower fluid cans- 
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Figure 7. Infrared images of the flow pattern for different air heights and different viscosities of 
silfcone oil. (a) and (b) used a 91 cS silicone oil. (c) and (d) used a 142 cS silicone oil. (a) and 
(c) had a 3 mm air height, (h) and (d) had a 20 mm air height. 


ing interfacial-tension-driveii flow in the low^er fluid. To explore this, a set of four 

experiments was performed. m i-o- ^ • 

In all of the experiments, a liquid aspect ratio of 2.0 was studied. Two different air 
heights (3 and 20 mm) and two different viscosities were used. When the air height of 
3 mm was used, the flow patterns did not change with viscosity but the temperature 
difference across the liquid did increase proportionally, indicating that convection 
was controlled liy the liquid phase. The liquid convection pattern also agreed with 

calculations and the experimental result is seen pictorially in figures 7a, c. 

When a deeper air height was used, the temperature difference across the liquid at 
the onset of convection did not change substantially between the experiments that 
employed different fluid viscosities. This indicated that convection was controlled by 
the dynamics in the air layer. It may be argued that air convection, when dominant, 
acts like buoyancy-driven convection between two rigid conducting plates as the 
lower liquid is much more viscous and more conductive than the air above it. A 
comparison was therefore made between the measured temperature drops across 
the air for the deeper air heights and the numerical calculations of Hardin et at. 
(1990). The experimental and theoretical results compared remarkably well and the 
flow pattern predicted theoretically also compared favourably with the experimental 
results. This confirmed our hypothesis that deep air heights interact with the lower 
liquid and drive thermally coupled flow through the interface. 

This phenomenon of thermal coupling may not be as important m LEVB because 
the crucible is often heated from above. However, this may be much more applicable 
to other important processes, such as evaporation and drying of films. 


3. Future work 

The research of convection in multiple fluid layers has revealed and continues to 
reveal many interesting phenomena. However, further work is needed to elucidate 
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some of the details more fully in a realistic system. One of the more interesting areas 
involves analysing some of the basic instability phenomena in bounded containers. 
To do this, two-fluid-layer numerical models that take into account realistic no-slip 
conditions will be necessary. With a proper model some of the effects, such as the 
Rayleigh-Taylor instability, the Gershuni-Zhukhovitskii instability and oscillations 
between thermal and viscous coupling, can be studied for containers with small 
aspect ratios. 

To date few experiments have been performed in small aspect-ratio containers. 
As was demonstrated in the codimension-two point experiments, new and interest- 
ing dynamics are present in small aspect-ratio containers, which are not present 
in large aspect-ratio containers. It would be interesting to show the interaction o 
codimension-two points with such instabilities as the oscillations between thermal 
and viscous coupling. Additionally, several of the phenomena discovered with theo- 
retical models have yet to be shown in experiments. Two examples are the Gershuni 
Zhukhovitskii instability and the oscillations shown by Rednikov et al. (1998). 

By investigating some of the basic physics of multilayer convection, we obtain 
a better understanding and appreciation for the liquid-encapsulated crystal-growth 
process and other fluid materials processing problems where temperature gradients 
are employed. Further research into this field should lead to improvements m such 
an important industrial process. 
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Discussion 

J. R. Helliwell {Department of Chemistry, University of Manchester, UK). Ref- 
erence has been made to oscillatory convection flow patterns in fluids and that this 
is known to cause defects (dislocations and fault lines) in crystal growth. My own 
particular research interests include the growth of protein crystals for X-ray crystal- 
structure analysis and how the quality of protein crystals can be improved, and 
thereby exploited, for higher resolution X-ray crystallographic data collection. I have 
been using CCD and interferometry diagnostic monitoring of protein crystal growth, 
and have seen benefits of microgravity if the crystals do not move, and if the mother 
liquor is not subject to convection (including Marangoni convection). The benefits 
of these conditions manifest as reduced crystal mosaicity and likewise a lack of, or 
only a few, mosaic blocks in X-ray topographs of crystals in such cases. In his exper- 
iment. how can Dr Johnson be sure that it is specifically oscillatory flow patterns 
that especially caused defects in his type of crystal? 

D. Johnson. We cannot be sure that oscillatory convection is always responsible 
for defects in crystals. However, research cited by Hurle (1994) has indicated that 
oscillatory behaviour generated through double diffusion is the cause of striations 
along the growth axis in directional solidification. The point of this paper, however, 
is to show that oscillatory behaviour need not arise merely from opposing forces that 
are seen in thermo-soiutal, otherwise known as double-diffusive, convection. Such 
oscillatory behaviour can arise by geometrical effects. Indeed as the crystal grows, 
the aspect ratio of the liquid phase changes and there are certain aspect ratios where 
the energy states may coexist leading to codimension-two points that can cause 
oscillatory convection. 

S. K. WJlson {Department of Mathematics, University of Strathclyde, Glasgow, 
UK). I complement the authors on a penetrating investigation of a complicated 
physical situation. As I understand it. they have found examples of slow oscillations 
between two different steady flow patterns in the vicinity of codimension-two points 
calculated theoretically using linear theory for the onset of steady convection in a 
finite-sized container. May I ask if truly oscillatory (rather than quasi-steady) con- 
vection is ever observed, and if it would be possible to undertake the same kind of 
investigation for the onset of oscillatory convection? 
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D. Johnson. Yes. under certain circumstances, we believe that oscillatory convection 
can be observed in liquid-gas bilayer experiments. Theoretical calculations that were 
done indicate the absence of such convection at the onset. In that ca.se what is the 
origin of oscillations in our experiments? The answer lies in the fact that a Hopf 
bifurcation lies in the vicinity of the onset but only in the post-onset region. The Hopf 
bifurcation mode or oscillatory mode was excited by the presence of codimension- 
two points. These points were generated by the fact that at certain aspect ratios 
two competing flow states coexist and in a manner of speaking the system wants to 
choose between the flow states leading to continual oscillations. 
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In this paper our most recent research results on natural convection in a closed cylinder, where our 
interest focuses on pattern structure dependence on aspect ratio and on temperature-dependent 
viscosity, are summarized. The main results are (a) the experiments on the onset pattern and 
conditions for pure Rayleigh convection in circular cylinders compare favorably with linearized 
stability results of Hardin et ai [Int. J. Num. Methods Fluids 10, 79 (1990)], as well as 
three-dimensional nonlinear calculations made by us; and (b) experiments and nonlinear 
calculations indicate a variation of the patterns at and near the codimension two points when large 
temperature differences are introduced, so as to cause a substantial change in viscosity. © 7995 
American Institute of Physics. 


I. INTRODUCTION 

The physics of Rayleigh convection is well understood 
and the classical theory has been reviewed in the treatise by 
Chandrasekhar’ and in the recent book by Koschmieder." In 
this problem, a layer of fluid is heated from below, and one 
of the objectives is to determine the conditions for the onset 
of flow from an erstwhile quiescent state, as well as the 
associated pattern. The critical temperature difference for the 
onset of flow is given in terms of the Rayleigh number. This 
group expresses the ratio of the buoyancy force, which pre- 
cipitates the convection, to the effects that dissipate the flow, 
viz thermal and momentum diffusivities. In the theory, the 
critical Rayleigh number is seen as a bifurcation point from 
the unstable quiescent solution. The experimental determina- 
tion of critical Rayleigh numbers is often affected by the 
presence of vertical sidewalls and the review by Azouni^ 
deals with the modulation of sidewalls on convection. 

Stork and Miiller*^ performed experiments to determine 
the critical conditions at the onset of convection, and they 
compared their results with calculations of Charlson and 
Sani.^ Some of these calculations were recently corrected 
and extended by Hardin et al^ An important observation that 
can be made from these calculations is that at critical condi- 
tions two azimuthal modes can coexist at certain aspect ra- 
tios. These are called codimension two points, and it turns 
out that they are spaced well apart at small aspect ratios 
(radius/height), but occur more frequently as the aspect ratio 
increases. It is also seen from these calculations that as the 
aspect ratio increases the various bifurcation points corre- 
sponding to various azimuthal modes lie close to each other. 

The experimental verification of the calculations is the 
main reason for our interest in this problem. Most experi- 
ments were performed by applying a series of temperature 
differences across a bounded liquid layer and visualizing the 
flow to determine the onset state. It can be seen from the 
definition of the Rayleigh number that when experiments are 
conducted in deep liquid layers, the critical temperature dif- 
ferences become very small (scaling inversely with the cube 
of depth). However, the time constant to reach steady condi- 
tions increases quadratically. Consequently, we may ask 


whether any adverse effects can occur by lowering the liquid 
depth so as to get better temperature control across the liquid 
layer, thereby reducing the time constant so that frequent 
changes of experimental conditions can then become pos- 
sible. The main consideration is that the thermophysical 
properties can now vary considerably over the liquid height. 
While it is evident that the critical Rayleigh number will now 
deviate from the predictions of classical theory that assumes 
constant properties, it is not so clear how the patterns will 
change and in what regions of aspect ratio these changes will 
be prominent. The purpose of this paper is to provide the 
evidence, both numerically and experimentally, which illus- 
trates that viscosity variation does affect the position of the 
codimension two points. A part of this study is devoted to the 
numerical modeling of nonlinear three-dimensional (3-D) 
convection showing the patterns near the onset for fluids of 
constant thermophysical properties. Additionally, the effect 
of temperature-dependent viscosity on flow patterns is shown 
numerically along with experimental verification. 

II. THE EXPERIMENTAL APPARATUS 

A schematic of the experimental test unit is shown in 
Fig. 1 . Experiments were conducted with Dow Coming sili- 
cone oil sandwiched between two plates. The lower bound- 
ary was made of anodized aluminum (thickness=3 mm). 


Cooling Water in 



no. I. Schematic of experimental apparatus for Rayleigh convection. 
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TABLE I. The properties of the Dow Coming silicone oil used in the ex- 
penments. 


T "bottom ‘^bottom ’^^top) Z / H 


Thermal conductivity 
Specific heat 
Density 

Thermal diffusivity 
Kinematic viscosity 
at 35 “C 

Thermal expansion 
coefficient 


0.001 588 W/cmK 
1 .463 J/g K 
0.968 g/cm^ 
0.001 cm Vs 

0.697 cmVs 

0.000 96/K 


while the upper boundary was an optically transparent sap- 
phire (thickness = 10 mm). Flow visualization was made pos- 
sible by aluminum flakes that were dispersed throughout the 
oil. The sidewalls were precisely machined Lucite. The oil 
was derated under vacuum so as to remove most of the en- 
trapped air. It was then loaded into the test section, taking 
care to avoid the presence of bubbles. The dimensions of the 
test section were measured with a micrometer and were ac- 
curate within 0. 1 mm. The lower side of the anodized alumi- 
num plate was in contact with a uniformly stirred water bath, 
which, in turn, was heated from below by an electric resis- 
tance heater. The bath was filled with water and bubbles 
were removed. It was felt that a uniformly stirred bath under 
the test section was better than direct contact of the alumi- 
num plate with a resistance heater, as we were interested in 
avoiding externally imposed thermal signatures on the test 
section. The upper sapphire plate was in contact with a con- 
tinuously flowing water bath. The temperature difference 
across the test section was assumed to be the temperature 
drop between the lower surface of the aluminum plate and 
the upper surface of sapphire, as both aluminum and sap- 



initial condition: velocity = 0 
temperature = 3D disturbance 


calculate velocity 


^ satisfy 

continuity no- 

scquationV^ 


calculate temperature 


same ^ 
as last 
Jtcration?> 


correct 

pressure 



FIG. 3. The 3-D temperature disturbance used in the calculations. 


phire have large conductivities compared to the silicone oil. 
The thermophysical properties of the Dow Coming silicone 
oil used are given in Table 1. Lucite has same thermal diffu- 
sivity as silicone oil within a percent and is 12%-15% more 
thermally conductive. Thus, the temperature gradients in the 
sidewalls and working fluid, i.e., silicone oil are very close, 
leading to a “conducting sidewall’* thermal boundary condi- 
tion, We will observe in a later section that the main results 
of this paper do not depend significantly whether the side- 
wall thermal condition is insulating or conducting. It may be 
noted that the thermal diffusivities of aluminum (0.97 cm^/s) 
and highly conductive sapphire (0.1633 cm^/s) are large 
compared to silicone oil, and so the time constants for a 
temperature change to establish steady conditions in these 
solid plates were quite short. The temperature differences 
were measured with a pair of resistance temperature detec- 
tors attached to the water bath side of aluminum and sap- 
phire plate. They were calibrated by an Omega 700 ther- 
mistor with the accuracy for temperature difference 
measuring within a 0.05 °C. The dynamic viscosity of sili- 
cone oil versus temperature is given as 
fjL—O.QM 471 ^^^g/ cm s. This came from the viscos- 

ity measurement at several temperature levels. The dynamic 
viscosity was measured in a coaxial cylinder, normal open 
cup Haake RV-12 and the viscosity measurements were ac- 
curate within 4%, assuming Newtonian behavior. The water 
temperature in the bath above the sapphire plate was con- 
trolled at a fixed point, while the water temperature in the 
bath below the aluminum plate was changed according to the 
heating patterns for each experimental run. The temperature 
differences across the test section were controlled with a PID 


TABLE II. The onset Rayleigh number and flow pattern indices reported by 
Hardin et ai^ 


FIG. 2. The flow chart of the calculation procedure. 


RIH 

Critical Ra 

Flow pattern 

1 

2260 

(0,2) 

1.8 

1835 

(L3) 

2.5 

1781 

(0,4) 

0.75 

2592 

(LI) 
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(a) R/H=l 




(c) R/H=2.5 


(d) R/H=0.75 


FIG. 4. The flow patterns from the calculations with constant viscosity: (a) R!H — \, Ra— 2270; (b) /?///— 1.8, Ra— 1900; (c) RIH 2.5, Ra 1800, and (d) 
/?///=0.75, Ra=2600. 


controller, which was tuned for the experimental conditions. 
The heating pattern in all of the experimental runs consisted 
of several segments at different values of temperature differ- 
ences around the onset point, and each of them had a slow 
ramping-up period and a constant temperature difference pe- 
riod greater than four times the thermal time constants of the 
test section. These time constants were calculated from the 
thermal diffusivity and the largest length scale in the test 
section. 


ill. THE MODEL 

The calculation was performed in a cylindrical coordi- 
nate system where z is the vertical coordinate. The governing 
equations are of the well-known Boussinesq form, except 
that the viscosity is taken to be a function of temperature. On 
scaling the governing equations in a manner similar to Har- 
din et the following definition of Rayleigh number (Ra) 
is arrived at Ra=/3g AT H^Ikv, Here, H is the depth of the 
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OOP 

(b)R/H=i.8, (1,3) 



IC )C )( )( ) 

(c) R/H=2.5, (0,4) 




FIG. 5. The schematics of flow patterns from Hardin et al.:^ (a) R/H= 1, (b) R/H = 1.8. (c) RIH=2.5, and (d) RIH=0J5. Here X indicates the falling fluid 
while O indicates the rising fluid. 


fluid and A 7 is the temperature drop across it, g is the gravi- 
tational constant, P is the thermal expansion coefficient, 
which is positive, and v and k are kinematic viscosity and 
thermal diffusivity at a reference temperature. 

The density was taken to be linearly dependent upon 
temperature, based on information from the Dow Coming 
Company that supplied the oil. The other thermophysical 
properties were obtained from the Dow Coming Company 
and are reported to be insensitive to temperature changes 
within the range used in this study. It is for this reason that 
only the viscosity variation with temperature was considered 
important in this study. 

The boundary conditions assumed no slip and no flow on 
velocities at the rigid surfaces, with vanishing temperature 
perturbations at the horizontal boundaries. Along the vertical 
boundaries most of the calculations assumed insulating sides, 
even though thermally conducting sidewall conditions that 
used a constant temperature gradient in the solid vertical 
walls would be more appropriate. As mentioned later, the 
conclusions of this investigation did not change very much 
when insulating vertical sidewall conditions were used. 

The SIMPLE^ (Semi-Implicit Method for Pressure Linked 
Equations) algorithm was employed to carry out the nonlin- 


ear calculations. The flow chart in Fig. 2 shows the proce- 
dure of the calculation. 


IV. DISCUSSION OF THE RESULTS 

In order to test the accuracy of the model, the computa- 
tions were performed for the classical Rayleigh problem, by 
which is meant a closed, rigid, and circular cylinder contain- 
ing a “Boussinesq’’ fluid heated from below and insulated 
along the vertical sides. As the model is a nonlinear 3-D 
simulation, no attempt was made to obtain the exact bifurca- 
tion point. Instead, the Rayleigh number was an input param- 
eter to the model. A 3-D disturbance on the linear distribu- 
tion of temperature was introduced to generate the onset of 
convection. Figure 3 shows the temperature distribution as 
an initial condition for the calculations. The initial tempera- 
ture distribution was horizontally uniform at a value that was 
linear in the z direction, except the two areas marked, where 
it was equal to the bottom temperature. The resulting steady 
pattern was calculated from this disturbance. When the Ray- 
leigh number was close to the bifurcation point reported by 
Hardin et ai!" the predicted pattern was identical to what 
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FIG. 6. The comparison of flow patterns between experiments and the calculations with a variation ot viscosity; (a) R!H 1.5, AT 19.0 C, (b) RfH 1.8 
Ar^ 18.5 X; and (c) R/H = 2.5, \T= 1 8.0 °C. 


was reported by them. A quantitative comparison of the ve- 
locity and temperature field was not possible, as Hardin 
et aL^ did not report this aspect of their study. The results ot 
our calculations are given in Figs. 4(a)~4(d) and those of 
Hardin et al.^ are given in Table fl for a variety of aspect 
ratios (/?///), where R is the radius of the test section. Note 
that two graphs, with the top view and cross section, are 
presented for each case in Fig. 4. The location where the 
cross-section graph was made is indicated by a dashed line in 
the top view graph. The input Rayleigh numbers for each 
case are shown in the caption of Fig. 4. Table II contains the 
bifurcation point and the pair of indices that shows the flow 
structure at onset. The first index represents the azimuthal 
mode number while the second index represents the maxi- 
mum number of roll cells across the diameter of the cylin- 
drical test section. For example, a flow mode (0,2) represents 
the schematic shown in Fig. 5(a), while mode (1,3) is seen in 
Fig. 5(b), mode (0,4) in Fig. 5(c), and mode (1,1) in Fig. 
5(d). Again, two graphs are presented for each case in Fig. 5. 
In the top view graph, the circles indicate the flow coming up 
from the plane of the paper, while the crosses show the flow 
going down into the plane. It may be observed that the re- 


sults of the calculations as shown in Fig. 4 are in good agree- 
ment with those given in Table II and Fig. 5. 

The next phase of the study involved the experimental 
determination of the patterns and comparison with the nu- 
merical modeling. Experiments were conducted in a variety 
of aspect ratios, and since the viscosity varied across the 
depth of the layer, our interest was focused on the patterns 
that would evolve as a result of this complication. Calcula- 
tions were done so as to incorporate the viscosity variation. 
The results of this comparison are shown in Figs. 6(a)-6(c) 
for aspect ratios of 1.5, 1.8, and 2.5, respectively. The calcu- 
lation results are shown in two graphs tor each case: (1) a 
3-D flow pattern at the outer surfaces of a part of the cylin- 
der, and (2) a flow plan at a constant “c" plane near the top. 
Experimental results are presented by the photographs taken 
from the top of the test section. In these photographs, the 
white regions are the locations where mainly horizontal 
flows occurred, while the black regions were occupied by 
mainly vertical flows. In each case, the comparison is excel- 
lent. The experimental runs were repeatable and the same 
working fluid was used in each experimental run. The flow 
visualization was made possible by the introduction of alu- 
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FIG. 7. The comparison of computational results and experimental results at the aspect ratio near the codimension two point; (a) /?//y=1.6, H = 1.2 mm. 
Ar=3,9 °C; (b) R/H= 1.6 //=4.2 mm, AT= 19.9 X. 


minum flakes, which were in the submicron range size into 
the liquid so as to form a uniform slurry. The difference in 
the numerical results shown in Fig. 4 and Fig. 6 is that the 
latter assumed a viscosity variation. Even though the actual 
experiment simulated the conducting sidewall case, the 
model reflected the insulating sidewall case. In fact, a sepa- 
rate calculation with the conducting sidewall case for aspect 
ratios 1.5 and 1.8 showed very little discernible change in 
patterns from the result, where an insulating wall condition 
was used. The comparison given in Fig. 6 is very encourag- 
ing, and there seems to be very little to distinguish these 
results from the perfect Rayleigh problem results of Figs. 
4(b) and 4(c). Because the viscosity variation was taken into 
account, the comparison between experiments and model are 
shown for the actual temperature difference (AT) used in our 
experiment. 

This raises the question of what effect, if any, the vis- 
cosity variation has on the patterns at onset. It is surmised 
that the location of a codimension two point must be predict- 
ably affected. A codimension two point is one where two 
modes coexist for the same critical Rayleigh number. In 
other words, it is the aspect ratio at which a switch in pat- 
terns takes place from one state to another. Straughan^ has 
shown in the case of an unbounded layer of fluid that the 
critical wave number is reduced when a viscosity variation 
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with temperature is assumed. This observation should there- 
fore immediately translate into a shifting of the codimension 
two point to the right, i.e., toward larger aspect ratio, at least 
for the case of a cylinder with artificial conditions of vortic- 
ity free sidewalls. It is to be noted that the unfolding of the 
critical Rayleigh number versus wave number curve into the 
relationship between Rayleigh number and aspect ratio has 
been used in the past by Charlson and Sani^ and Rosenblat,'^ 
among many others. To test the hypothesis that the codimen- 
sion two point would be shifted to the right if a fluid of 
varying viscosity were used, a number of experiments were 
conducted in various aspect ratios with tall and short heights 
of liquid layers. The corresponding calculations were made 
with the assumption of insulating sidewalls as the position of 
the codimension tw'o point depends only slightly on this 
change. For example, Hardin et ciL^ report that the codimen- 
sion two point for conducting walls occurs at an aspect ratio 
of 1.59 and at an aspect ratio of about 1.58 for the case of 
insulating sidewalls. According to the calculations of Hardin 
et a flow mode of (0,2) should occur immediately to the 
left of the computed codimension two point, while a mode of 
(1,3) .should occur immediately to the right. Figure 7 shows 
the comparison of patterns between modeling and experi- 
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ment for an aspect ratio of 1.6 for a tall fluid height of 7.2 
mm, with a temperature difference of 3.9 °C and for a short 
fluid height of 4.2 mm, with a temperature difference of 
19.9 The comparison between modeling and experiments 
is again very good and a remarkable change in patterns was 
observed. The tall fluid set into a mode (1,3) and the short 
fluid set into a mode of (0,2), indicating that the codimension 
two point had shifted to the right when a larger temperature 
difference was introduced, such that the variation of viscos- 
ity became large enough to make the shifting distinguishable. 
This observation thereby agreed qualitatively with our fore- 
cast, which, in turn, was based on the calculations of 
Straughan^ for a laterally unbounded fluid. 

V. SUMMARY 

Comparison of experimental results with results of de- 
tailed 3-D calculations for Rayleigh convection in a closed 
cylinder was made. Very good agreement was obtained be- 
tween theory and experiment. The viscosity variation did not 
affect the onset flow patterns significantly from the constant 
viscosity case, except at the vicinity of the codimension two 
points. Based on the theoretical prediction of Straughan® for 
a laterally unbounded fluid, it was conjectured that the codi- 
mension two points would shift toward larger aspect ratios if 
viscosity variation were to be considered. This conjecture 
was shown to be correct in our study and once again we 
obtained excellent agreement between experiment and the 
corresponding numerical modeling. 
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